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PREFACE 


This document contains the proceedings of the Workshop on Computational Methods for 
Crashworthiness held at NASA Langley Research Center, September 2-3, 1992. The workshop 
was jointly sponsored by the University of Virginia Center for Computational Structures 
Technology and NASA. Workshop attendees came from government agencies, the aerospace and 
automotive industries, energy laboratories, and universities. The objectives of the workshop were 
to assess the state-of-technology in the numerical simulation of crash and to provide guidelines for 
focused future research leading to an enhanced capability for numerical crash simulation. 

Certain materials and products are identified in this publication in order to specify 
adequately the materials and products that were investigated in the research effort In no case does 
such identification imply recommendation or endorsement of products by NASA, nor does it 
imply that the materials and products are the only ones or the best ones available for the purpose. 
In many cases equivalent materials and products are available and would probably produce 
equivalent results. ■ 


Ahmed K. Noor 

University of Virginia Center for Computational Structures Technology 
Hampton, Virginia 


Huey D. Carden 

NASA Langley Research Center 

Hampton, Virginia 
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INTRODUCTION 

The numerical simulation of the nonlinear response of aircraft structures and other 
components subjected to large impact loads has been the focus of intense efforts because of the 
pressing needs for collision protection, or crashworthiness, of transportation vehicles (aircraft, 
cars, trains and ships), and nuclear reactors and containment vessels. Several software systems 
are currently available for impact analysis of structures (e.g., KRASH, DYNA, MSC/DYTRAN, 
NIKE, DYCAST, SUPERWHAMS, PAMCRASH, PRONTO and McALOG/RADIOSS). 

Crashworthiness is already a design requirement for many vehicles and is likely to become 
more important with increasing demands for safety considerations in structural design. Major 
advances are needed in our analysis capabilities for predicting the structural response of vehicles 
subjected to impact loads and the subsequent response of human occupants. To this end, there are 
a number of technology needs and related tasks that must be addressed by the research community 
to enhance the state-of-the-art in computational methods for crashworthiness. 

The joint NASA/University of Virginia workshop held at NASA Langley Research Center, 
September 2-3, 1992 focused on the status of computational methods for crashworthiness and the 
current and future needs for further development of this technology. The list of the pacing items 
given in this introduction was compiled from a number of participants. Research efforts to address 
these needs can guide the design of future transportation vehicles in the following two ways: 1) by 
providing better understanding of the phenomena associated with crash, thereby identifying the 
desirable energy-absorbing design attributes; and 2) by verifying and certifying crashworthy 
designs, and making low-cost modifications during the design process. 

COMPUTATIONAL NEEDS FOR CRASHWORTHINESS 

The technology needs identified by the participants can be grouped into the following five 
major headings: 1) understanding the physical phenomena associated with crash; 2) high-fidelity 
modeling of the vehicle and the occupant during crash; 3) efficient computational strategies; 4) test 
methods, measurement techniques and scaling laws; and 5) validation of numerical simulations. 
For each of the aforementioned items, attempts should be made to exploit the major characteristics 
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of high-performance computing technologies, as well as the future computing environment. The 
five primary technology needs and related tasks are described subsequently. 

1. Understanding the Physical Phenomena Associated with Crash 

This includes understanding a) the mechanics of large dynamic deformations of structures, 
including the effects of frictional contact; b) the effects of inertia forces and of material strain-rate 
sensitivity on the dynamic response; and c) damage initiation and progression during crash. For 
the occupant, the factors that can be correlated with the level of injury or death (e.g., dynamic 
response index, head injury criteria, force in the lumbar spinal region) need to be identified. The 
modeling details required to capture the different phenomena associated with the structural 
response during crash need to be identified. 

2. Hi gh-Fidelity Modeling of Vehicle and Occupant 

The reliability of the predictions of the response of the structure and occupant during crash 
is critically dependent on: a) the accurate characterization and modeling of material behavior, b) 
high-fidelity modeling of the critical details of the vehicle and occupant (e.g., seat, fasteners, and 
the human anatomy); and c) modeling of the frictional contact between the vehicle and the 
impactor, as well as between the different parts of the vehicle, including the need for accurate 
material constitutive models and properties for foam, padding, and textile composites, especially 
the strain rate sensitivities, for modeling of the seat/occupant interaction. 

3. Efficient Computational Strategies 

The effective use of numerical simulations for predicting the vehicle response during crash 
requires strategies for treating phenomena occurring at disparate spatial and time scales, using 
reasonable computer resources. The strategies are to be based on using hierarchical (multiple) 
mathematical models in different regions of the vehicle to take advantage of the efficiencies gained 
by matching the model to the expected response in each region. To achieve the full potential of 
hierarchical modeling there should be minimum reliance on a priori assumptions about the 
response. This is accomplished by adding adaptivity to the strategy. The key tasks of the research 
in this area are the following: 
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a) Rational selection of a set of nested mathematical models for different regions of the 
vehicle and discretization techniques for use in conjunction with the mathematical models. The 
mathematical models should include damage models. 

b) Efficient numerical algorithms for simulating frictional contact and other local 
phenomena, such as stiffness reduction due to damage and redistribution of load paths in damaged 
components. 

c) Automated (or semiautomated) coupling of different mathematical/discrete models. 

d) Sensitivity analysis to assess the sensitivity of crash response to each of the material and 
geometric parameters used in the computational model. 

e) Criteria for adaptive refinement (or derefinement) of the mathematical and discrete 
models. 

f) Stable and efficient iterative procedures and numerical algorithms for use in conjunction 
with adaptive model refinement. 

4. Test Methods. Measuring Techniques and Scaling Laws 

The effective coupling of numerical simulations with experiments requires a high degree of 
interaction between the computational analysts and the experimentalists. This is done at three 
different levels, namely: 1) laboratory tests on small specimens to obtain material data; 2) 
component tests to verify computational models and to determine empirical structural properties 
which can be used in hybrid experimental/numerical models; and 3) full-scale (or scale model) tests 
to validate the computational model and assess the need for model improvements. 

New test methods and measurement techniques are needed to study progressive failure, as 
well as soil and water impact. The influence of specimen size or scale factor on structural response 
is not well understood. Thus, testing of geometrically similar sub-scale models is not possible, 
until the scaling laws governing the phenomenon are understood. In particular, scaling laws are 
needed which account for the material behavior including elastic properties, failure initiation and 
ultimate strength, structural and topological details, as well as the loading characteristics. 
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5. Validation of Numerical Simulations 


In addition to validating the numerical simulations by component and full-scale tests, a 
number of carefully selected benchmark tests are needed for assessing new computational 
strategies and numerical algorithms. These benchmark tests would provide a measure of 
confidence in new codes, or added functional capabilities to existing codes. They could also serve 
as a basis of code comparisons for efficiency and accuracy for modeling of impact problems 
involving large structural deformations in short time durations. 

RELATED TASKS 

In addition to the aforementioned technology needs, some related tasks need to be 
addressed before numerical crash simulations can have a significant impact on the design process. 
These are: 1) development of software for automated (or semiautomated) model (and mesh) 
generation; 2) pre- and post-processing software for efficient input and reduction of numerical 
simulation data; 3) use of advanced visualization technology; and 4) adaptation of AI tools 
(knowledge- based/expert systems and neural networks) to crash simulation systems. 
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Highlights of the Workshop 


Ahmed K. Noor 

Center for Computational Structures Technology 
University of Virginia 




Objectives and Format 
of Workshop 


The numerical simulation of the nonlinear response of structures subjected to large impact loads has 
been the focus of intense efforts in recent years. Despite these efforts, major advances are needed in 
various areas of the technology before numerical simulations can be routinely used to satisfy 
crashworthiness design requirements of transportation vehicles. The objectives of the present workshop 
(Fig. 1): to assess the state-of-technology of computational methods for crashworthiness; and to identify 
current and future needs for further development of the technology. 

The workshop includes presentations and two panels. The presentations are included in the 
proceedings to illuminate some of the diverse issues, and to provide fresh ideas for future research and 
development. 


Objectives 

• To assess the state-of-technology in the numerical 
simulation of crash 

• To identify future directions of research 


Format 

• Presentations 

• Panels 

• Panel 1 - Computational Needs for the Accurate 
Simulation of Crash. Moderator: Ed Fasanella 

• Panel 2 - Experimental Needs and the Coupling 
Between Experiments and Numerical Simulatons. 
Moderator: Huey Carden 

• Proceedings 


Figure 1 
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Assessment of the 
State-of-Technology 


The first aspect of assessing the state-of-technology is to assess our understanding of the physical 
phenomena associated with crash. Some of the issues that affect these physical phenomena are listed in 
Fig. 2. These are: mechanics of large, dynamic deformations; material-level damage mechanisms, 
damage growth and subsequent structural failure; material strain-rate sensitivity; contact/impact conditions 
and friction modeling; energy absorbing characteristics of structures; and the influence of specimen size or 
scale factor on structural response. 


Understanding of Physical Phenomena 
Associated with Crash 

• Mechanics of large dynamic deformations 

• Damage and failure mechanisms 

• Effect of material strain-rate sensitivity 

• Effect of inertia forces 

• Contact/impact conditions 

• Friction modeling 

• Energy absorbing characteristics of structures 

• Scaling laws 


Figure 2 
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Assessment of the State-of-Technology 

(Cont'd.) 

Current Capabilities 


The second aspect of assessment of technology is that of current capabilities for numerical 
simulation of crash (Fig. 3). These capabilities include computational material models and damage 
mechanisms, modeling of structural details such as joints, seats, fasteners and occupant; efficiency of 
currently-used computational strategies for handling spatial discretization and temporal integration, 
frictional contact/impact conditions, hierarchical, global-local and adaptive refinement facilities; and 
assessment of software systems currently used for impact analysis of structures. 


• Computational Material Models 

Constitutive models and material data, friction models, 
damage mechanisms and failure models 

• Level of Details 

Modeling of the structure, seat, fasteners and occupant 

• Efficiency of Computational Strategies 

- Spatial discretization and temporal integration 

- Hierarchical, global-local, multilevel ana adaptive strategies - 
interfaces between models 

- Contact/impact algorithms 

• Current Software Systems 

KRASH, MSC/DYTRAN, SUPER WHAMS, WRECKER, 
DYCAST, DYNA, NIKE, PRONTO, PAM CRASH, 
McALOG/RADIOSS 


Figure 3 
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Future Directions of Research 


Three factors should be taken into account in identifying future directions for research (Fig. 4): 

1) characteristics of future transportation vehicles and their implications on design requirements for 
crashworthiness; 

2) future computing environment and computing paradigm; and 

3) recent and future developments in other fields of computational technology, which can be adapted 
to numerical simulation of crash. 

Two of the important research tasks are: 

1) validation of numerical simulations and selection of benchmark tests for assessing new 
computational strategies and numerical algorithms. The standardized tests would provide a measure of 
confidence in added functional capabilities to existing codes, or in new codes; and 

2) treatment of uncertainties in material properties, geometry, boundary conditions, and operational 
environment through probabilistic analysis, stochastic modeling and sensitivity analysis. 


• Characteristics of future vehicles and their implications on 
design requirements for crashworthiness 

• Impact of emerging and future computing environment (high- 
performance computers, multimedia workstations, advanced 
visualization technology) 

• Impact of developments in other fields of computational 
technology (e.g., CFD, computational mathematics) 

• Validation of numerical simulations and effective coupling 
with experiments (Benchmarks) 

• Treatment of uncertainties in material properties, geometry, 
boundary conditions, spatial and temporal distribution of 
loading and operational environment (probabilistic analysis, 
stochastic modeling and sensitivity analysis) 

Figure 4 
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ADAPTIVE METHODS FOR NONLINEAR STRUCTURAL 
DYNAMICS AND CRASHWORTHINESS ANALYSIS 


Ted Belytschko 
Northwestern University 
Evanston, Illinois 


ABSTRACT 

The objective of this talk is to describe three research thrusts in crashworthiness analysis: 

1) adaptivity 

2) mixed time integration, or subcycling, in which different timesteps are used for different parts of 
the mesh in explicit methods 

3) methods for contact-impact which are highly vectorizable. 

The techniques are being developed to improve the accuracy of calculations, ease-of-use of 
crashworthiness programs and the speed of calculations. The latter is still of importance because 
crashworthiness calculations are often made with models of 20,000 to 50,000 elements using explicit time 
integration and require on the order of 20 to 100 hours on current supercomputers. 

The methodologies will be briefly reviewed and then some example calculations employing these 
methods will be described. The methods are also of value to other nonlinear transient computations. 
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OUTLINE 


• Adaptive mesh procedures in nonlinear 
analysis: why, how, and what is the 
status 

• Subcycling (mixed time integration) 

• New highly vectorizable methods for 
contact impact which are well suited to 
adaptive methods 

Figure 1 


PREDICTION 

The 1990’s will be the decade of adaptivity, 
adaptive mesh refinement 
adaptive targeting 
adaptive organization objectives 

Figure 2 
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PREDICTION 


There are three types of adaptivity, which are known by the letters r, h, and p. These letters are 
mnemonic letters and refer to how the refinement is achieved. In r methods the nodes are relocated. In h 
methods, refinement is achieved by reducing the element size h. In p methods, refinement is achieved by 
increasing the order p of the element interpolance. 


TYPE OF MESH ADAPTIVITY 


r — method 



relocate nodes 


h — method 


adapt element size h 


p — method 



adapt order p of element interpolants 


Figure 3 
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ADAPTIVITY IN NONLINEAR FEM 


Adaptive methods are particularly useful in nonlinear problems such as crashworthiness because 
nonlinear response is often characterized by localization. In the areas of localized response more 
refinement is needed. When standard method is used, the user of the program must refine the mesh where 
he anticipates this localized deformation. Therefore, different meshes must be developed for different 
loadings. For example, in car crash, different meshes must be developed for frontal and rear impact, side 
impact, and overturning. This can be quite expensive from the viewpoint of manpower. 


Why are adaptive methods particularly important in nonlinear 
problems? 

Modes of failure of structures 

i. buckling, particularly with formation of hingelines 

ii. localization 

iii. fracture 


All of these involve local phenomena whose location cannot be 
determined at the outset of a simulation. 


Figure 4 
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COMMENTS ON ADAPTIVITY FOR SHELL AND 
CRASHWORTHINESS PROBLEMS 


In comparing the different types of adaptivity for nonlinear structural dynamics problems such as 
crashworthiness, the following advantages, which are marked by a plus sign (+), and disadvantages 
which are marked by a minus sign (-), can be attributed to the various types of methods. From this study 
we concluded that the h-method was the most suitable method for adaptivity in crashworthiness. 


r — method 

- large elements cannot represent shape of shell 

+ most accuracy with given NDOF compared to h 

- history diffusion 

- elements become distorted - decreases accuracy 
+ easiest data structure 

p — method 

- awkward in nonlinear dynamics; no good lumped mass 
+ easy data structure 

h — method 

+ relatively effective 
+ no distortion of elements 

- moderately complex data structure 


Figure 5 
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TYPES OF ERROR INDICATORS 


Error indicators are an important ingredient in adaptive methods since they are to a guide the refinement 
of the mesh. Error indicators are classified by Oden in the following classes: residual, interpolation, and 
post-processing. In the work we are doing, we are using projection error criterions, a post-processing 
type, because they are very easy to implement and are quite effective for low-order elements. 


1. Residual: Compute residual in governing equations and use its 
norm or use it to drive an element or local enriched solution. 

a) Explicit: Evaluate a norm of the residual. 

b) Implicit: Use residual to drive a local or element error 
equation. 


2. Interpolation Methods: Estimate magnitude of derivatives of 
higher-order than contained in finite element space. 


3 . Projection (postprocessing) Methods: Obtain a smoothed solution 
and compare to finite element solution; sometimes called L2 
projection methods. 


Figure 6 
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ADAPTIVE SCHEMES FOR TRANSIENT AND 
NONLINEAR PROBLEMS 


based on constant resource approach 

1 . Advance the solution n time steps 

2. Compute element error indicators 0 e 

3. Sort 0e 

4. Fission elements with 0 e > tol fusion 
Fuse elements with 0 e < tol fu sion 

5. Repeat the last n time steps with new mesh (optional) 

6. go to 1 


Note: If n is too small or tol fiss *o n too close to tol fusi( >n j W e 
encounter "churning" which degrades accuracy. Our recent 
experience shows 5 is quite important. 


Figure 7 



REMARKS ON H-ADAPTIVITY 


Constraints (or slave nodes in explicit methods) must be introduced 
at nodes where a large element has two or more neighbors on one 
side to enforce compatibility; easy in vector methods, awkward in 
matrix methods. 


Usually a group of contiguous elements should be fissioned 
simultaneously because fissioning a single element does not provide 
much enrichment; only one new free node. 


In wave propagation problems, change in element size can cause 
spurious reflections. 


Usually mesh gradation is limited to 1 -irregular meshes: large 
element cannot have more than 2 small neighbors on any side; see 
Devloo, Oden and Strouboulis (1987). 


Data structure with fission and fusion is complex, particularly for 
real engineering meshes; see Belytschko, Wong and Plaskacz, 
Computers and Structures, 33(4-5), 1989, pp. 1307-1323. 


Figure 8 
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MIXED TIME INTEGRATION 


In h-adaptive meshes, large variety of element sizes are found. When explicit methods are used of 
such meshes, the tunestep is reduced dramatically by the presence of small elements. Therefore methods 
called mixed time integration (or subcycling) are being used. 


Motivation : in explicit integration with same At over 
entire mesh, stiffest element sets At. also called subcycling, 
explicit-explicit partitions; 


example 


h_ 

10 





[?2 




II hK- 

LN 


— i 

h k- 
/■ 


B 


Atcrit = min (-) 


c = wave speed 


for A 



h 

c 


for AuB 


At 


h 

crit = 10c 


so AuB is lOx as expensive as A 


Figure 9 
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Mixed Time Integration 


Integrate each element or subdomain with At cr it using an interface 
treatment that preserves stability + consistency. 


In example 


2 4 


' 2 

3 

4 

5 


9 

1 0 


1 3 


to — H I h H 


integrate element 1 and nodes 1 to 4 with 



elements 2 to 10 and remaining nodes with 

At-| 

cost savings: ~ 90% 

In adaptive methods, large range of stable time steps is unavoidable, 
so subcycling is crucial for efficiency. 


Figure 10 
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CONTACT-IMPACT 


The modeling of contact-impact is very important in the simulation of crashworthiness. However, 
contact-impact algorithms often require more than fifty percent of the running time of a crashworthiness 
code because they are not easily vectorized. Therefore we have developed a pinball algorithm which is far 
more highly vectorizable. 


Contact-impact is an important phenomenon in crash analysis, e.g., 

1. engine impact with body, fire wall 

2. wheel impact with inner fender 

3. contact of collapsing surfaces 

Most contact-impact algorithms require many different branches. 



The branch of the algorithm which is activitated depends on which 
surface is penetrated; there are special branches for edges, etc. 


Figure 1 1 
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PINBALL PENALTY ALGORITHM 


T. Belytschko and M. O. Neal, International Journal for Numerical 
Methods in Engineering, 31, 1991, pp. 547-572. 


Interpenetration and inteij>enetration rate g are 
computed on pinballs inserted in elements. 



Enforces contact-impact conditions on spheres embedded 
in elements. 

As h — » 0, impenetrability is enforced. 

Algorithm is simple and highly vectorizable. 


Figure 12 
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Salient Features of Algorithm 


Radius of pinball is determined by equivoluminal expression 




471 


Pinballs are classed by body; for single-surface slideline, smaller R 
needed. 


Interpenetration has occurred when 


d. < R + R 
y i j 

g = d.. 

& u 


Pinball forces are equally transferred to all nodes of associated 
element (a surface node option available). 


The pinball method automatically places pinballs on outside 
elements by using assembled surface normal algorithm. 


Figure 13 
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EXAMPLES OF NONLINEAR ADAPTIVE 
COMPUTATIONS 


Nonlinear, transient computations with an explicit nonlinear finite 
element program WHAMS using h-adaptivity and pinball for contact 
impact; see Belytschko and Yeh (1992). 


An L2 projection on the strain invariants was used to calculate an 
error estimate. 


A commercial version of this program is available from: 

KBS2, Inc. 

455 Frontage Road 
Burr Ridge, IL 60521 
(708) 850-9444 
Fax (708) 850-9455 


Figure 14 



TWO-LEVEL ADAPTIVE MESH OF CYLINDRICAL PANEL 


shows f ? n h- adaptive solution of a cylindrical panel which is impulsively loaded Notice that the 
and hingeTn« “ thC S “ PPO,t ’ Where lhere is Severe plasdc b '" di "* ^formation, 



Two-level adaptive mesh of cylindrical panel. 

Figure 15 
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This figure shows a comparison between solutions obtained by h-adaptivity and those obtained using a 
very fine mesh and a coarse mesh. As can be seen, the adaptive solution compares well to the fine mesh 
solution. The differences in the displacements obtained by the coarse mesh and the fine mesh are not 
large, but for some of the stresses and strains, significant improvement is obtained by the use ot 
adaptivity. 



Time 


Displacement of point A 
in cylindrical panel ( MAXLEV=2) 



Tint 


Strain Exx of element B 
in cylindrical panel (MAXLEV=2) 



Tent 



Strain Eyy of element B Stress Sxx of element B 

in cylindrical panel (MAXLEV=2). in cylindrical panel (MAXLEV=2) 



Stress Syy of element B in cylindrical panel (MAXLEV=2). 


Figure 16 
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This shows results for an S-beam which is impulsively loaded. Again, significant differences occur in 
some of the strains for a coarse mesh solution as compared to an adaptive or fine mesh solution. 



T - shape cross section 

Material number 2 ( Table 6 ) 


Geometry of T-shape cross section beam. 



Axial displacement Dz of point A 
in T-beam (MAXLEV=1). 



Strain Exx of element B 
in T-beam fMAXLEV=l). 



Axial velocity Vz of point A 
in T-beam (MAXLEV=1) 



Strain Eyy of element 
in T-beam (MAXLEV=1) 


Figure 17 
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This shows the evolution of the mesh for the S-beam; note that the h-refmement occurs at a comer 
where local buckling takes place. 




Time = 1.28 ms ; 673 clem 



Time=1.92 ms ; 688 clem 


One-levei adaptive mesh of T-beam. 


Figure 18 
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This is the same problem with a higher level of adaptivity. Far more elements are placed in the region 
of local buckling. 



Time=0. ms; 1 147 elem 



Time=0,32 ms ; 2794 clem 



Time = 0.96 ms ; 2581 elem Time=1.6 ms ; 2212 elem 


Two-level adaptive mesh of T-beam. 
Figure 19 
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This again compares displacement in strains for coarse mesh, fine mesh, and adaptive solutions. 
Again the adaptive solutions agree very closely with the fine mesh solution. 



Axial displacement Dx of point A 
in T-beam (MAXLEV=2). 



Strain Exx of element B 
in T-beam (MAXLEV=2) 



in T-beam (MAXLEV=2). 



Axial velocity Vz of point A 
in T-beam (MAXLEV=2) 



Strain Eyy of element B 
in T-beam (MAXLEV=2) 



Stress Syy of element B 
in T-beam (MAXLEV=2). 


Figure 20 
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This is a solution of a box beam which has an initial velocity as shown with an attached mass at the 
back. This problem is often considered a model for crash analysis. The solutions for fine mesh, coarse 
mesh and adaptive meshes are shown; the adaptive solution agrees well with the fine mesh. 



Geometry: L = 0.1500 m Initial condition: V= 15.64 m/sec 

a = 0.0300 m Material number 3 (Table 6) 

t = 0.0015 m 


Box beam problem. 



Velocity of point A in the box-beam (MAXLEV=1). 


- 1.000 10 * 

-1.500 10* 

-1000 10 * 

■1500 10* 

-3.000 10* 

-3.500 10* 

•4.000 10* 

4.500 10* 

6 a 001 0.002 0.003 0.004 0.005 a 006 
Tus 

Reaction force of box-beam (MAXLEV=1). 
Figure 21 
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This shows a timing for a full car model which is shown on the next page. It is solved with full 
contact-impact and subcycling. The important thing to notice is that subcycling gives a speedup of 1.7 and 
that the effective element cycle time on a CRAY-YMP here is 12 microseconds. 


Timing 

FULL-CAR MODEL 


Elements: 

17,297 

Mass (kg): 

1,880 

Time steps: 

78,274 

80 msec simulation 


CRAY-YMP 


Without subcycling: 
128 elements/block 

7.63 hrs 

Element cycle time: 

20 jisec 

With subcycling: 
64 elements/block 

4.39 hrs 

Effective element cycle time: 

12 jisec 

Speedup: 

1.7 
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Figure 22 



wlnBB - von-mlses shell material - subcycle 
time ■ 0.000E+00 



Figure 23 
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winBB — von— mises shell material — subcycle 

time ■ 0. 000E+00 



Figure 24 


32 


uiin8B - von-mises shell material - subcycle 
time ■ 4.001E+01 



Figure 25 
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Implementation and Stability 

Time steps are assigned to nodes and element blocks automatically. 
Elements are sorted by At£ nt 


At crit < At crit < < At crit 

e 2 n 


Elements are arranged in blocks so that time steps of adjacent 
elements have integer ratios.* 


Blocking of elements is necessary to take advantage of vectorization. 


For analysis of stability, see Belytschko and Lu, ASME publication 
edited by G. Hulbert, et al, 1992. 


*A new algorithm which does not require integer ratios has recently 
been developed (Belytschko and Neal, Computer Methods in 
Applied Mechanics and Engineering , 31, 1989, pp. 547-570). 


Figure 26 
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REMARKS AND CONCLUSIONS 


• H-adaptivity is a promising technique for 
simulating nonlinear structural response and 
structural failure. 

• Improves accuracy. 

• Simplifies model preparation. 


• Subcycling and advanced contact-impact methods 
such as the pinball method can improve efficiency 
of explicit dynamic codes and is essential with h- 
adaptivity. 


• Improved error criteria are needed for adaptive 
methods for nonlinear solid mechanics. 


Figure 27 
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INTRODUCTION 


The Landing and Impact Dynamics Branch of NASA Langley Research Center has been involved in 
impact dynamics research (Fig. 1) since the early 1970's. For the first ten years, the emphasis of the 
research was on metal aircraft structures in both the General Aviation Crash Dynamics Program (Refs. 
1-13) and the Controlled Impact Demonstration (CID) Program, a transport aircraft program culminating 
in the controlled crash test of a Boeing 720 aircraft in 1984 (Refs. 14-16). Subsequent to the transport 
work, the emphasis has been on composite structures with efforts directed at understanding the behavior, 
responses, failure mechanisms, and general loads associated with the composite material systems under 
crash type loadings. Considerable work has been conducted to address the energy absorption 
characteristics (Refs. 17-20) and it indicates that composites can absorb as much if not considerably 
more energy than comparable aluminum structures. However, due to their brittle nature, attention must 
be given to proper geometry and designs to take advantage of the good energy absorbing properties 
while providing desired structural integrity. Achieving the desired new designs often requires an 
understanding of how more conventional designs behave under crash type loadings. 

The purpose of this paper is to present a review of the composite impact dynamics research being 
conducted at NASA Langley Research Center. Examples are presented of experimental and analytical 
data to illustrate the activities in the four program elements of the composite research. 



Figure 1 
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IMPACT DYNAMICS RESEARCH FACILITY 


The research was conducted by personnel at the Langley Impact Dynamics Research Facility 
(Fig. 2) using test equipment located at the installation. The Impact Dynamics Research Facility 
(IDRF), originally the Lunar Landing Facility used by the astronauts during the Apollo Program for 
simulation of lunar landings, has been modified to allow crash tests of full-scale aircraft under controlled 
conditions. The aircraft are swung by cables from an A-frame structure which is approximately 400 ft 
long and 230 ft high. The impact runway can be modified to simulate other ground crash environments, 
such as packed dirt, to meet a specific test requirement. 

Each aircraft is suspended by cables from two pivot points 217 ft above the ground and allowed 
to swing pendulum- style into the ground. The swing cables are separated from the aircraft by 
pyrotechnics just prior to impact. The length of the swing cables determines the aircraft impact angle 
(from 0 degrees (level) to approximately 60 degrees). Impact velocities can be obtained up to 65 mph 
(governed by the pullback height). Variations of aircraft pitch, roll, and yaw can be varied by changes in 
the aircraft suspension harness attached to the swing cables. Data from onboard instrumentation are 
transmitted through an umbilical cable hard wired to the control room at the base of the A-frame. 
Photographic data are obtained by onboard, ground-mounted, and A-frame mounted cameras. 

Maximum allowable weight of the aircraft is 30,000 lb. Reference 21 provides complete details of the 
facility and test techniques for full-scale aircraft testing. 



Figure 2 
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COMPOSITE IMPACT DYNAMICS RESEARCH PROGRAM ELEMENTS 


The program elements of the Composite Impact Dynamics Research Program are illustrated in Fig. 
3. Currently, efforts in crash dynamics research are in four areas: (1) development of a data base to 
understand the behavior, responses, failure mechanisms, and general loads associated with both 
conventional and innovative concepts using composite material systems under crash type loadings; (2) 
analytical studies/development relative to composite structures; (3) studies in scaling of composite 
structures under static and dynamic loads; and (4) full-scale tests of metal and composite structures to 
verify performance of structural concepts. 

The overall goal of the research efforts is to gain a fundamental understanding of composite crash 
behavior and to formulate improved structural designs to meet performance, integrity, and energy 
absorption requirements. Examples of experimental testing relative to each program element will be 
highlighted in the paper. Analytical examples associated with the program elements will be presented at 
the appropriate time within the discussions of the experimental tests rather than in separate sections. 



Figure 3 
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TYPICAL TRANSPORT FUSELAGE STRUCTURE 


If one examines a typical transport fuselage structure as shown in Fig. 4, it becomes readily 
apparent that frames are one of the most important components used in the construction. As a 
consequence, the initial efforts in our studies were evaluation of individual frame components. The 
approach of studying simple structural elements and then moving to combinations of these elements in 
more complex substructures has been taken in the development of a data base on the dynamic response 
and behavior of composite aircraft structures. With this building block approach, more complex 
subfloor structures fabricated from the simpler components for both static and dynamic testing are 
discussed later in the paper. 


Floor bum Itypl 


Seat track (typ) 
Floor panel 



Pressure bulkhead 

Win) center box structure 


Fairings 


Section no. 4 
wing center section 
Keelson 


Section no. 2 


• Section no. 1 
flight station 

s Forward pressure bulkhead 


Figure 4 
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COMPOSITE FUSELAGE FRAME CONCEPTS 


Various cross-sectional shapes for fuselage frames are used in metal aircraft and are often 
proposed for composite structures. Figure 5 presents sketches illustrating four of the more common 
geometries, I-, J-, C- and Z-cross sectional shapes. Several circular frames using these shapes were 
fabricated for testing. To add out-of-plane stability to the frames (with the exception of the Z-section 
frames), 3-1/2 inch wide skin material was added enhancing the ease of testing of both symmetrical and 
other nonsymmetrical sections. The skin, a [±45/0/90] 2s lay-up sixteen ply (.08 inches) thick, was 
cocured with the 6 foot diameter frames. The frames were constructed in two different heights, 1-1/2 
inches and 3/4 inches, to investigate the effect of frame height on behavior and responses. 


C-section 

Z-section 



Figure 5 
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DYNAMIC LOADING BEHAVIOR OF COMPOSITE Z-FRAME 


One of the first geometries to be studied under static and dynamic loadings was the Z-cross 
section. The six-foot diameter frames were constructed of 280-5HA/3502, a five harness satin weave 
graphite fabric. The height of the frame was 3 inches with a total width of 2.25 inches and about 0.08 
inches thick. Lay-up of the frames was quasi- isotropic. Initial tests were with 360 degree frames made 
from four 90 degree segments joined with splice plates. Additional tests were conducted with half 
frames since the top half of the complete frames were undamaged in the tests. 

Figure 6 presents results from the dynamic studies of the response of the Z-frames under 
approximately a 100 Ibm floor loading at an impact velocity of 20 fps. As noted in the figure, the splice 
plates joining the segments of the frame are 45 degrees up die circumference from the point of impact. 
Also, complete failures (fractures) of the Z-section frames occurred at the bottom and approximately 60 
degrees from the bottom. Potentially, it appears that the presence of the splice plates may have 
influenced the locations by moving the top failure points up a few degrees to about the 60 degree 
locations. 


Failure locations 


20 fps 20 fps 
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COMPOSITE FUSELAGE FRAME TEST 


As a result of earlier tests with the Z cross-section circular frame, a 3.5 inch wide skin material 
was added to the I-, J-, and C- frame concepts to increase the torsional stiffness of the cross-sections and 
limit out-of-plane rotations and deformations. The skin, a [±45/0/90]2s lay-up sixteen ply (.08 inches) 
thick, was cocured with the 6 foot diameter frames. Data in Fig. 7 are for a lay-up of the frame of 
[±45/45/90/03 ] s . Both the skin and the frame were fabricated with AS4/5208 graphite-epoxy material. 

Figure 7 shows a typical set-up of a composite fuselage I-frame in a 120 000- lbf loading 
machine prior to a quasi-static test. TTie purpose of the study was to evaluate the effect of floor 
placement on the structural response and strength of the circular fuselage frames constructed of 
graphite-epoxy composite material. A steel I-beam was attached horizontally across the composite 
frame at the diameter position to simulate the floor. The horizontal floor positions were designated by 
the included angle measured between the ends of the floor attachments about the center of curvature of 
the frame. For example, the frame with the floor at the diameter is designated the 180° floor since the 
arc is 180° between the attachment points. 

A vertical compressive load was applied to the composite fuselage frame through the simulated 
floor beam and the lower platen of the load machine. Special clamps (see lower right of Fig. 7) were 
used to bolt the I-beam to the composite frame and a 170° F melting point metal was poured into the 
small gap between the clamp and the frame to eliminate possible motion in the joint. As shown at the 
bottom left and top right of the figure, additional tests were conducted where the floor location was 
moved to produce 120° and 90° arcs. In each test the specimen was loaded at a rate of 500 lbf/minute 
up to a maximum of 1000 lbf. 



Figure 7 
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FINITE ELEMENT MODEL OF I-FRAME 


To gain an understanding of the physics of behavior, the experimental research of structures 
under crash loadings is generally accompanied by analytical prediction or correlation studies whenever 
feasible. Various finite element codes with capabilities for handling dynamic, large displacement, non- 
linear response problems of metal and composite structures were used as tools in the research efforts. 
The analytical results presented in this paper were generated with a nonlinear finite element computer 
code called DYCAST (DYnamic Crash Analysis of STructures (Ref. 22) developed by Grumman 
Aerospace Corporation with principal support from NASA and the FAA. 

Using the beam element from the DYCAST element library, the I-beam model illustrated in Fig. 
8 was formulated. The combination of outer skin and I-frame were modeled with only I-beam 
elements. Since the skin lay-up provided less stiffness than the lay-up of the 1-frame, the skin width 
was reduced by the ratio of the computed stiffnesses of the skin to the I-frame. As a result, the 3.5 inch 
skin width was reduced to approximately the same 2.5 inch width as the bottom flange of the I-frame 
itself. Thus, the resulting model consisted of straight ISEC elements with the bottom flange and skin 
combined to be 0.16 inches in thickness with only the material properties of the I-frame being used in 
the model. Symmetry was utilized and thirty-nine I-beam elements were used for the beam model. 

Both a force and a moment loading was applied to the top end of the frame. The static analytical load 
was increased linearly to a maximum load of 1000 lbf in 50 pound increments. The analytical results of 
the model are compared to the experimental data in the following section. 


APPLIED 
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Figure 8 




EXPERIMENTAL ^ND ANALYTICAL STRAIN COMPARISON FOR STATIC 
LOADING OF THE COMPOSITE I-FRAME 


the 0° or ground SiSStodSn Tt^wo s3d^ P m°f * ^ mdude: (a) maximum drains were * 
between ± 55°-60° from the bottom contact area- 3t symmetric loc ations 

die same "sea gull" shape as measured in the ext)eriment- P lS F Skin strain dlstn bution exhibits 

distributions were noted for the inner flange of the frame’ as rvl ^^Hverted circumferential strain 

in the magnitudes of the analyticKS f?5* d “ the experiment. The agreement 

considered excellent. expenmental strains and the shape of the distributions are 

composite frame were: (a) to^tenhe mag^md^of the^af 0 ?" changblg the floor position in the 
"sea gull" shape of the distribution under vertical loadbi^rS ^° me ^. but not the common, general 
strain distribution to occur in the frame segment below the floor 22!?" 111 generaI " sea gul1 ” sha P ed 
increase rhe effective globa, s„ srS&^Tn^'^S ^° ased . 


Strain 


180° Floor 



Circumferential location, deg 


120° Floor 



Circumferential location, deg 


Strain 


90° Floor 

4 



Circumferential location, deg 


° Experimental skin 
D Experimental flange 
— Analytical skin 
■** Analytical flange 


Figure 9 
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COMPOSITE SUBFLOORS WITH AND WITHOUT SKIN 


As indicated previously, the approach of 

combinations of these elements id aircraft structures. The approach parallels 
data base on the dynamic response andbe aircraft programs. Consequently, three composite 

^ “ S d ' SCUSSMi ab0VC ' 

Figure 10 is a photograph of the skeleton iXted* 

of the single Z-section frames similar to bondine methods. Aluminum floor beams tied the 

the three frames through metal clips and sec ?^ d ^ . , f the su bfloor. Notches in the frames allowed 

^ ^ &L tests were confer, with the 

subfloors. 



Figure 10 
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COMPOSITE SUBFLOOR BEHAVIOR-SPECIMEN WITHOUT SKIN 


For the three composite subfloor specimens used for impact studies, two static and two dynamic 
tests were conducted on the subfloors. With the skeleton subfloor, a static and a dynamic test to 
destruction whs conducted. ^Vith the skinned subfloor, a non-destructive static test followed by a 
dynamic test to failure was conducted. 

Figure 1 1 shows the locations of fractures of the skeleton subfloor after an impact test onto a 
concrete surface at 20 feet per second. In the dynamic test of the skeleton subfloor, fractures were 
produced at notches in the frames. The locations, shown in Fig. 1 1, were also near the point of impact 
(about 11 degrees because of the splice plate) and at two other locations up the circumference of the 
frames (45 degrees and 78 degrees) and involved all three frames for a total of 15 fractures. The impact 
energy exceeded the energy absorbed by the local fractures and the floor bottomed out in the impact. 


Skeleton subfloor 

Failure locations 
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COMPOSITE SUBFLOOR BEHAVIOR-SPECIMEN WITH SKIN 


Figure 12 presents impact results for the skinned subfloor after an impact of 20 feet per second. 
Points of failure of the frames in this specimen are indicated in the figure. Again the points of failure are 
at/near the impact point (within 12 degrees) and circumferentially at about 56 degrees up both sides of 
the frame on the middle and back frame and 45, 12 and 22.5 degrees on the front frame. It was observed 
that the subfloor impacted first on the front area which possibly explains the 12 and 22.5 degree 
fractures being different from the other locations. Again all three frames were involved in the failures. 
Some delamination of the frames from the skin was evident but the skin remained intact. 


Skinned subfloor 
Failure locations 
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COMPARISON OF FRAME BEHAVIOR WITH SUBFLOORS 


The determination of the effect of the floor location on the structural response of fuselage frames 
has aided in the understanding and prediction of full-scale subfloor or fuselage response to crash 
loading. For example. Fig. 13 shows a comparison of the normalized experimental dynamic strain 
distribution on the flange of the skeleton subfloor and the skin location corresponding to the flange 
position of the skinned composite subfloor specimens with the analytical I-frame strain. The results 
from the simple frame show a strong similarity to the response of the more complex subfloor structures. 
The structures share in common the generally circular or cylindrical shape, the vertical loading 
situations, and under vertical loads have strain (moment) distributions which have maximums at the 
point of loading and at approximately ±45° to ±60°, depending on boundary conditions, around the 
circumference from the ground contact point. Analytical results show the same distribution with 
maximums corresponding to the experimental locations. Failures of the subfloor structures were noted 
between these same 45° to 60° circumferential locations in the dynamic tests (see Ref. 23). 


Composite Structures 


Normalized 

experimental 

strain 



Figure 13 
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METAL TRANSPORT FUSELAGE FAILURE BEHAVIOR 


To relate previous metal transport failure behavior to the current studies and observations with 
composite structures, data from tests of metal aircraft sections to support transport aircraft research 
efforts are included. NASA Langley Research Center conducted drop tests of two 12-foot long fuselage 
sections cut from an out-of-service Boeing 707 transport aircraft to measure structural, seat and occupant 
responses to vertical crash loads, and to provide data for nonlinear finite element modeling. One section 
was from a location forward of the aircraft wings and one was from aft of the wing location (Refs. 24 
and 25). The sections were loaded with seats, anthropomorphic dummies, data acquisition system pallet, 
power pallet, and camera batteries to test not only structural, seat, and occupant responses but also to test 
equipment to be used in the full-scale transport crash conducted later. 

Structural damage locations of the transport aircraft structures resulting from the 20 fps drop tests 
are shown in Fig. 14. The damage to the transport sections was confined to the lower fuselage below the 
floor level. Under the vertical impact of 20 fps, all of the frames ruptured near the bottom impact point. 
Plastic hinges formed in each frame along both sides of the fuselage at about 50 degrees up the 
circumference from the bottom contact point. The upward movement of the lower fuselage was 
approximately 22-23 inches at the forward end and 18-19 inches at the rear for the section taken from 
forward of the wing location, whereas in the section from aft of the wing location the crushing was about 
14 inches forward and 18 inches in the rear. Although the aircraft structures are metal and the failures 
discussed above involve plastic deformations with some tearing of the metal rather than brittle fractures, 
the general observed failure pattern and locations for the transport fuselage sections are noted to be quite 
similar to the results of the composite frames and subfloors discussed herein. 


Metal transport sections 
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SUMMARY OF OBSERVED FAILURE BEHAVIOR 


The response behavior determined during the studies of full-scale aircraft sections, fuselage 
frames, and subfloors are summarized in Fig. 15. The figure shows normalized moment distribution on 
a representative frame of the various specimens and the failure locations which were noted from static or 
dynamic tests. The visual impression is quite striking among the various specimens. It is suggested that 
from the results of simpler frames to the more complex subfloors and full-scale sections, a strong 
similarity is evident in the failure behavior of the structures. The structures share in common the 
generally circular or cylindrical shape, the normal loading situations, and what appears to be a similar 
pattern of failure behavior. Analytical models of frame structures under vertical loads have moment 
distributions which have maximums at the point of loading and at approximately 45 to 60 degrees 
(depending on boundary conditions) around the circumference from the ground contact point. Failures 
of the structures were noted at these same locations. Such observations can help dynamists gain a better 
understanding of what to expect from such structures in crash loading situations, can guide designers of 
new structures to better account for the vertical crash loads, and allow better energy absorption to be 
included in the new designs. Additionally, the observations can help analysts better model the aircraft 
structures for predicting the failure responses and behavior under crash situations. 
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FULL-SCALE COMPOSITE AIRCRAFT CONCEPT 


As was illustrated in Fig. 3, one of the four program elements of the Composite Impact 
Dynamics Research in the Landing and Impact Dynamics Branch of the Structural Dynamics Division at 
NASA Langley Research Center is full-scale testing of aircraft structures under crash loading conditions. 

Two full-scale composite general aviation aircraft structures, two complete wing sets, and 
landing gears have been obtained for testing. As shown in Fig. 16, the full-scale test aircraft is an 8 
place airplane with twin Pratt and Whitney engines (650 hp each) powering a pusher type propeller. 

The gross takeoff weight is approximately 7200 lbm with empty weight being approximately 4000 lbm. 
Overall length of the plane is 38 feet with a wing span of 39 feet. The structure of the test aircraft 
represents a composite skin/frame type fuselage construction concept with the exception of the interior 
floor structure which consists of aluminum beams on which the seat rails and seats are mounted. 

Design support testing is underway to replace the existing floor structure in one of the two aircraft with 
an energy absorbing concept constructed with composite materials. The retrofit approach is, of course, 
necessitated and as a consequence, special tests have been undertaken to develop the replacement floor 
concept to assure that structural behavior and failure loads and modes of failure are achieved in the 
concept prior to inclusion in the aircraft. 



Figure 16 
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ENERGY ABSORBING BEAM DESIGN FOR COMPOSITE AIRCRAFT SUBFLOOR 


A preliminary experimental program is being conducted to study the behavior of energy 
absorbing composite spars for aircraft subfloors. The study, which is a part of a wider full scde aircraft 
test program, examines the efficiency of replacement floor structures for an existing all composite 
fuselage aircraft shown in the lower right of Fig. 17. The efforts are a continuation of previous research 
(Refs. 7 and 13) dealing with crashworthy metal aircraft subfloor structures. A typical section of the 
composite fuselage with the original subfloor structure is shown in the center of the figure. As shown in 
the center figure, the four spars that support the seat rails are aluminum, whereas the rest of the subfloor 
structure is graphite composite. Static tests of such a subfloor section have shown that the existing 
structure is too stiff and too strong for cushioning loads resulting from crash speeds in the neighborhood 
of 30 fps, as recommended by the Part 23 of the Federal Aviation Regulations. Thus, the objective of 
this study is to design and test a retrofit subfloor structure that would provide the desired cushioning 
(less than 20 g of occupant load) at crush speeds of approximately 30 fps. In particular, the four 
aluminum sparsj are to be replaced by composite sine wave beams. The sine wave composite beam 
concept, shown at upper left, has been examined previously with encouraging results (Ref. 19). 


EA SPAR 



Figure 17 
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PROPOSED COMPOSITE SUBFLOOR STRUCTURE 


A schematic of the proposed subfloor structure is shown in Fig. 18. An ideal beam design should 
contain two flanges for improved stiffness which, in addition, would offer two crush initiators, one in 
juncture between the web and the flange. However, because of the double curvature of the fuselage and 
the retrofit nature of the problem, the beams under construction are limited to a single flange on the top 
which contains the single crush initiator. 



Figure 18 




STATIC RESPONSE OF INITIAL COMPOSITE SPAR DESIGN 


A preliminary spar beam with an inverted "L" cross section has been fabricated according to the 
suggested geometric details in Ref. 19. The stacking sequence and composite systems were selected 
such that the stiffness of the beam transverse to the longitudinal axis was half that of the aluminum 
spars. Both graphite and aramid reinforced epoxy have been used in the stacking sequence; 
(±30°gr./±45°kevlar fabric/0°2 gr.) s . A series of static and dynamic tests have been conducted to 
evaluate the overall performance of the spar design. A typical virgin and crushed spar section is shown 
in Fig. 19 (a) and a load/displacement plot of a quasi-statically loaded composite spar is shown in Fig. 

19 (b). Note that, while the section crushed progressively, thus absorbing a large amount of energy, the 
ultimate load and the sustained crushing loads of approximately 1200 and 650 lb/in respectively are far 
too high as opposed to 200 - 300 IbAn of desired load. 
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Figure 19 
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DYNAMIC RESPONSE OF INITIAL COMPOSITE SPAR DESIGN 


A series of dynamic tests were carried out on a 14 ft drop tower shown in Fig. 20 (a) to simulate 
a 30 fps mass drop equal to the corresponding weight of the seat and the occupant which would load the 
spar section. The mass for a 12" long spar section was 184 lbm. A typical load/time plot from a 
dynamic test is shown in Fig. 20 (b). Note that, while the value of the sustained crushing load was 
comparable to the one obtained in the corresponding static tests, the ultimate load was much lower. The 
dynamic results, shown in Fig. 20 (b), indicate that there was no dynamic rate effect up to the 30 fps 
impact velocity of the test as compared to the static data and that the loads from the dynamic test of the 
preliminary beam design is also much too high for a human occupant. 




(b) 

Figure 20 
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CONSIDERATIONS FOR NEW COMPOSITE SPAR DESIGN 


The previous results together with the experience gained in the fabrication of the first spars 
influenced the design of a new set of spars and tests. Some of the factors and constraints that controlled 
the design of the new spars included: 

(1) simplification of the fabrication process - eliminate unidirectional prepreg, 
complex angles, and hybridization 

(2) elimination, if possible, of the graphite reinforced material to improve ductility 
under dynamic loads- to ensure survivability of the flange under normal landing 
loads 

(3) improvement of the web bending stiffness to ensure adequate longitudinal spar 
bending stiffness following the loss of the flange 

(4) reduction of the ultimate and the sustained crushing loads to less than 300 lb/in 
(and greater than 200 lb/in) - to improve cushioning 

(5) symmetric loading - apply load symmetrically from the flange to the web to 
improve global stability of the spar- web and improve stroke efficiency 

It was found that more design goals could be met with a sandwich construction and that some of 
the additional complications associated with the fabrication of the sandwich spars could also be offset 
with the simplification of the skin lay-up. Thus, a "T" section illustrated at the bottom of Fig. 21 was 
chosen instead of the original "L" section. A number of sandwich spar sections are being fabricated with 
fabric kevlar webs and hybrid flanges. A full test matrix is shown above the sketch of the "T" section. 
Static and dynamic testing of these sections will commence after specimens have been instrumented. 
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FINITE ELEMENT MODELING OF THE COMPOSITE AIRCRAFT 


As part of the overall program, along with the full-scale testing of the composite aircraft concept, 
finite element predictions of the behavior and loads of the aircraft under crash conditions will be 
conducted. The computer program KRASH (Ref. 26), which is a three-dimensional, hybrid, finite 
element modeling technique, will be used to predict response and loads of the test aircraft under selected 
impact parameters using the KRASH finite element aircraft model (or variation thereof) shown in Fig. 
22. As depicted in the figure, KRASH represents the structure of the aircraft as a combination of 
masses, beams, rigid connections, and external springs. 



Figure 22 
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SCALING OF COMPOSITE STRUCTURES 


The final research area in Composite Impact Dynamics Research is the scaling of composite 
structures. One activity is discussed herein, but Refs. 27 to 31 deal with other studies. 

Figure 23 shows some typical results of a study to isolate the factors responsible for scale effects 
in the tensile strength of angle ply graphite/epoxy composite laminates. Two generic ±45° lay-ups were 
studied, one with blocked plies and one with distributed plies with stacking sequences containing 
between 8 and 32 plies. The 8-ply laminate consisted of off-axis plies arranged in a (+45°/-45°)2S 
sequence, and was denoted the baseline or model stacking sequence. A high modulus, high strength 
brittle graphite-epoxy system (AS4/3502) was used to fabricate six "scaled-up" laminates with the 
following stacking sequences: (+45°n/-45°n)2S (blocked plies), and (+45%45°)2nS (distributed plies), 
where n = 2, 3, and 4. Tensile coupon specimens having four scaled sizes were constructed including 
full scale size, 3/4, 2/4, and 1/4, corresponding to n equal 4, 3, 2, and 1, respectively. Angle ply 
laminates are commonly used for damage tolerance in the surface of composite laminates where the load 
bearing plies are shielded against impact and fatigue loads. It is therefore important to understand the 
effect of specimen thickness and stacking sequence on the stress-strain response, the ultimate strength, 
and the mode of failure for this class of laminates. 

Results in Fig. 23 indicate that for increasing specimen size: (1) strength decreased for blocked 
ply laminates; (2) strength increased for distributed ply laminates; and (3) strength of distributed 
laminates is greater than blocked laminates for a given specimen size. The significance of these findings, 
beyond the scaling of structures issues, is that ASTM standard tests for determination of the in-plane 
shear stiffness and strength are based on ±45° angle ply testing, although exact specifications for the 
laminate stacking sequence are not stated. Results of this research show that the values of strength can 
vary tremendously depending on whether the laminate stacking sequence contains blocked or distributed 
plies. Also, the size of the laminate, especially the number of plies, is important. Recommendations to 
improve the standard testing practices have been made to the ASTM so that a meaningful shear strength 
value, independent of specimen size, can be determined from tensile tests on ±45° laminates. 
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CONCLUDING REMARKS 


The Composite Impact Dynamics Research Program at NASA LaRC focuses on generating a 
database for understanding composite structural behavior under crash loads, examines conventional and 
innovative metal and composite structures for meeting performance, integrity, and energy absorption 
requirements, analyzes/enhances analysis tools for composite applications, studies scaling effects in 
composite structures under static and dynamics loads, and conducts full-scale structures to verify 
performance of structural concepts. 

Typical examples of research in each of the program elements were presented to illustrate the 
research effort. Experimental and analytical results were presented which showed the effect of floor 
placement on the structural response of circular fuselage frames constructed of graphite-epoxy 
composite material. The results from the simple frame showed a strong similarity to the response of 
more complex subfloor structures. The structures share in common the generally circular or cylindrical 
shape, the vertical loading situations, and under vertical loads have strain (moment) distributions which 
have maximums at the point of loading and at approximately ±45° to ±60°, depending on boundary 
conditions, around the circumference from the ground contact point. Analytical results show the same 
distribution with maximums corresponding to the experimental locations. Failures of the subfloor 
structures were noted between these same 45° to 60° circumferential locations in the dynamic tests. 

A design support test program to develop a composite energy absorbing floor structure to replace 
metal floor in a composite aircraft concept was outlined and preliminary results presented. A 
preliminary spar beam with an inverted ”L" cross-section has been fabricated. A series of static and 
dynamic tests have been conducted to evaluate the overall performance of the spar design. A typical 
spar section crushed progressively absorbing a large amount of energy, but the ultimate load and the 
sustained crushing loads were far too high for human survivability. In the dynamic tests, the value of 
the sustained crushing load was found to be comparable to the one obtained in the corresponding static 
tests; however, the loads from dynamic tests of the preliminary beam design were also much too high 
for a human occupant. New design goals were established which should be met with sandwich type spar 
construction. 

Scaling results were presented which have had wide-spread influence on standard test practices 
for material properties. ASTM standard tests for determination of the in-plane shear stiffness and 
strength are based on ±45° angle ply testing, although exact specifications for the laminate stacking 
sequence are not stated. Results of this research have shown that the values of strength can vary 
tremendously depending on whether the laminate stacking sequence contains blocked or distributed 
plies. In addition, size of the laminate, especially the number of plies, is important. Recommendations 
to improve the standard testing practices have been made to the ASTM so that a meaningful shear 
strength value, independent of specimen size, can be determined from tensile tests on ±45° laminates. 

The Composite Impact Dynamics Research Program will contribute to the technology necessary 
for the development of improved composite structural aircraft concepts for energy absorption and 
enhanced passenger protection under crash loads. 
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STRUCTURAL CRASHWORTHINESS 


The key point of this definition is protection of the occupant during a crash event. The structure 
must be designed to absorb kinetic energy while controlling dynamic forces to within acceptable human 
tolerances and the maintenance of livable space around the occupant. 


• Definition: 

The ability of a vehicle to reduce dynamic forces experienced by 
occupants to acceptable levels while maintaining a survivable envelope 
around them during a specified crash event. 



STRUCTURAL CRASHWORTHINESS (Continued) 


There are three key points shown here. The first is that structural crashworthiness, made 
mandatory by government regulations, is now a design criteria for automobiles and helicopters. Because 
of this it is essential to assess vehicle crashworthiness during the design cycle. This involves simulation 
of nonlinear behavior of complex structures during impact which is computationally and theoretically 
intensive, involving state-of-the-art developments in computational mechanics. These developments are 
on-going. 


• Structural crashworthiness has become a design criteria for occupant 
carrying vehicles - especially for automobiles and helicopters. 

• A requirement exists to assess vehicle crashworthiness during the 
design cycle. 

• Numerical simulation of the dynamic response of vehicles subjected 
to impact loading is computationally intensive and has/will involve 
state-of-the-developments in computational mechanics. 
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ESSENTIAL REQUIREMENTS OF STRUCTURAL CRASHWORTHINESS 


The mechanics of structural deformation during a crash involves geometric nonlinearities due to 
large shape changes, material nonlinearity due to elastic-plastic behavior and nonlinearities due to variable 
contact/rebound of structural parts. 

Characterization of elastic-plastic behavior of metallic structures is reasonably well understood. 
Efficient computational algorithms have been developed and implemented into finite element programs. It 
should be emphasized that these algorithms are sensitive to actual material properties so that there is 
continuing need for experimental data. Impact of composite materials is a newer and challenging problem 
as the phenomena of failure differs from metals. Progressive failure and damage laws must be developed 
that describe fiber failure, matrix cracking, ply debonding and sublaminate buckling. These phenomena 
then must be implemented into simulation codes for crashworthiness of composite structures. Obviously, 
an extensive suite of tests must be performed ranging from coupon tests to subcomponents to full-scale 
section tests. 

Simulating variable contact between an impacting structure and external or internal surfaces is 
computationally intensive. Nevertheless, algorithms have been developed and implemented into a number 
of simulation codes. Contact friction between surfaces is usually treated using simple Coulomb friction. 


• Large elastic- plastic deformation with failure 

- accurate characterization of the constitutive material behavior 

- failure prediction for composite laminate construction 

• Variable contact/ re bound 

- contact with external surfaces 

- contact between internal parts 

- friction between contacted parts 
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ESSENTIAL REQUIREMENTS OF STRUCTURAL 
CRASH SIMULATION (Continued) 


• , I ?? ve ? prne [ ltsm com P uta tional algorithms are ongoing. Two areas come to mind that will increase 
the viability of crash simulation. These are temporal and mesh adaptivity, error estimates and parallel 
computations. The first area is necessary in order to assure that the discrete model is efficiently and 
accurately predicting structural impact behavior. The second is to provide the computing power to perform 
crash simulation on a routine and timely basis within the design cycle as well as to accommodate more 
detailed models that may be dictated by introducing adaptivity. 

In any occupant-carrying vehicle, there are parts designed specifically to absorb energy during a 
crash. In an automobile these can be rails that are designed to progressively crush. Very detailed models 
are required to simulate the accordian folds and internal contact of these parts although current full vehicle 
models include the rails in the entire model. The situation for aircraft and helicopters is equally complex, 
energy absorbing subfloor concepts have been developed that appear to be extremely difficult to model. 

I ne open question still remains. Is it possible to model these highly nonlinear regions in an accurate and 
cost-effective manner? Currently, crush data is developed from component testing and then implemented 
into a model as nonlinear springs. What about structures constructed of composite materials? 


• Accurate and efficient computational techniques 
? parallel computation 

? adaptive methods 

• Modeling capability for a variety of structural types 

- including special energy absorbing structural concepts 

- either metallic or composite 

- hybrid materials 


71 



MODELING STRATEGIES - PAST EXPERIENCE 


There are three distinct behavioral regions that must be considered when performing a crash 
simulation. 


• A model for crash simulation can be separated into three distinct 
modeling regions 

- linear 

- moderately nonlinear 

- extremely nonlinear 

• Pictorially this is shown on the next visual. 
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BEHAVIOR ZONE CHARACTERISTICS 


Yes, modeling the extremely nonlinear crushing regions is still a monster! 


BEHAVIOR ENCOUNTERED IN CRASH 
SIMULATION 



LINEAR MODERATELY EXTREMELY 

NONLINEAR NONLINEAR 
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BEHAVIOR ZONE CHARACTERISTICS - LINEAR 


Modeling the linear zone is obvious. Use as little computational resources as necessary. 


♦ Linear zone 

- elastic 

- small deflection 

• Modeled with 

- rigid bodies with lumped mass 

- relatively few elastic finite elements 

- substructure; most dot's omitted 
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BEHAVIOR ZONE CHARACTERISTICS - MODERATELY 

NONLINEAR 


These areas begin to be computationally expensive. Current technology in nonlinear structural 
mechanics is in place to adequately treat these regions. Allow for any possible global collapse modes. 


• Moderately nonlinear 

- elasto-plastic 

- large displacements on the global scale 

• Modeled with 

- nonlinear finite elements 

- allowance for possible global collapse modes 
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BEHAVIOR ZONE CHARACTERISTICS - EXTREMELY NONLINEAR 


These are the very difficult and computationally intensive areas to model. These areas often 
involve special energy absorbing components that exhibit large deformation on a local scale. The bottom 
line is that these parts require very detailed FEM models. 


• Extremely nonlinear 

- large deflection on a local scale, i.e., accordian folding of metallic 
structural components, crushing of subfloor structure in aircraft, 
local deformation and failure of composites 

- specially designed energy absorbing structure 

- crushable nonstructural parts 

- requires fine model (thousands of dof's) 
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BEHAVIOR ZONE CHARACTERISTICS - EXTREMELY 
NONLINEAR (Continued) 


In the past these areas were exclusively modeled with nonlinear spring elements. The explicit 
codes currently used make use of a detailed representation of highly nonlinear regions in automobile 
structural components. However, there still may be areas of a built-up structure that are currently not 
amenable to detailed modeling. This is particularly true for aircraft structures that involve complex energy 
absorbing floor concepts. Composites offer new challenges in detailed modeling of crushable 
components. More research is required in this area before viable crash simulation is feasible. 


• Modeled with: 

- Nonlinear spring elements 

Spring properties from test or other analysis require intimated 
understanding of the structural and material behavior 

- Detailed discrete representation?? 

Trade-off between "hybrid" model and detailed nonlinear finite 
element model 
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MODELING CRASHWORTHINESS OF COMPOSITE STRUCTURES 


In modeling the behavior of a composite airframe under impact conditions, the levels to which 
material failure must be included are far more detailed than for metallic airframes. In a traditional airframe, 
material failure modes which influence the structural failure process are primarily ductile crack growth and 
tearing and ductile rupture. Each of these is a distinct failure which renders a region of the structure 
incapable of transmitting any load. Just as important, because of the ductile nature of the material, 
significant impact energy can be absorbed by the material prior to failure. In laminates composed of 
graphite-epoxy lamina, failure at the material level is primarily brittle with very little prefailure energy 
absorption. For composite laminates, impact induced material damage modes such as matrix cracking and 
delamination lead to a structural material which still can transmit load, albeit at a reduced ultimate level and 
with a reduced stiffness and cross-sectional modulus. Thus, local material failures can have significant 
impact on the mode of structural failure. 

Because sublaminate material failure in composite laminates can propagate over large regions 
without leading to total laminate failure, the overall structural failure mode can be significantly altered. 
Buckling and collapse modes are highly dependent on the modulus and bending stiffness of the structural 
cross section. Zones of highly degraded material stiffness caused by fiber breaking/buckling, 
delamination, and matrix microcracking can thus cause structural failure modes significantly different than 
those seen in an undamaged structure of the same laminate construction. 


• Composite aircraft structures pose new computational and modeling 
challenges for crashworthiness design and analysis 

- brittle failure 

- diverse failure mechanisms 

- adhesive joints 

- three-dimensional effects in thick laminates 

• Failure modes in composite structures 

- buckling/collapse 

- fiber breakage (local) 

- delamination (structural) 

- microcracking (local and structural) 

• Material failure can change the buckling/collapse behavior 
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MODELING CRASHWORTHINESS OF COMPOSITE STRUCTURES 
CURRENT RESEARCH DIRECTIONS 


As a first step in modeling the crashworthiness of composite structures, a progressive failure 
algorithm was implemented into the DYCAST cede. This entailed introducing capability to describe 
classical laminate behavior, i.e., building a laminate from individual ply properties and orientation. This 
leads to the usual integrated material stiffness matrices; [A], the in-plane stiffness matrix, [D], the bending 
stiffness matrix, and [B], the matrix that couples the in-plane and bending behavior of the laminate for 
unsymmetric lay-ups. 

Initially three maximum strain failure criteria were implemented. These are described in the next 
figure. Currently, a strain based criterion based on the theory proposed by Christensen is being 
implemented and tested (see subsequent figures). 

Nonlinear material behavior was also investigated by implementing the Sandhu method for 
nonlinear composite materials. This method uses four uni-axial curves to define material behavior, two 
direct stress versus strain, shear stress versus shear strain and the variation of Poisson's ratio versus 
strain. Further investigation must be made of the effect of matrix material nonlinearities. 


* In response to the program directions set by the Impact and 
Dynamics Branch at NASA LRC, the following enhancements are 
being incorporated into DYCAST; 

• Ply-by-ply progressive failure laws 

- maximum strain 

- Christensen strain-based quasi three-dimensional 

- Tsai-Wu 

• Nonlinear material behavior 

- Sandhu method for nonlinear composite materials 

- nonlinear elastic shear response 
? matrix plasticity 
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MAXIMUM STRAIN FAILURE CRITERIA 


The table below describes the three maximum strain failure criteria initially implemented into the 
DYCAST code. The first two are for uni-directional composites and the third is for a fabric composite. 

The failure criterion chosen is checked at each finite element stress integration point in each ply of the built- 
up laminate. As indicated in the figure below, when fiber failure is detected, the moduli En and Poisson's 

ratio, Vi 2 as well as the stress. On, are set to zero. Similar procedures are followed for matrix and shear 
failure. Options 2 and 3 are similar to option 1 but add an induced coupling based on heuristic arguments 
that, for example, fiber failure induces shear ineffectiveness. Based on these considerations, the [A], [D], 
[B] matrices are reformulated and reflect the progressive softening of the damaged structure. 



Option 1 

(Unidirectional Composite) 

Option 2 

(Unidirectional Composite) 

Opt 

(Fabric Co 

Ion 3 
mposlte) 

Primary Failure Direction * 

Induced 

Additional 

Failure 

Zeroed Out 
Quantities 

Induced 

Additional 

Failure 

Zeroed Out 
Quantities 

Induced 

Additional 

Failure 

Zeroed Out 
Quantities 

Fiber Failure, c = p 

C 11 c u 

FAIL 

None 

E „’ V n 

Shear 

^ii ’ G\i j v , 2 
^11 ’^12 

Shear 

^12 

^ 11*^12 

Malrlx or Fiber Failure**, 

^22 ^22 

FAIL 

None 

F v 

L 22» M2 

o* 

Shear 

^22 <P\2 > V 12 
02 * ^12 

Shear 

F G v 

t, 2I* U l2» M2 
^22’ ^12 

Shear Failure**, 

7=7 

M2 M2 

FAIL 

None 

Gn 

T 12 

Matrix 

^12 »^22* ^12 
^12*^22 

Fiber 

Gil , Ell , E 22 , 

V 12 1 T 12 » Oil 

On 


* The subscript FAIL on &S denotes prescribed fail ure strain. Here, primary means the failure that Induces the additional failures. 
*# Note that with Option 2, primary matrix and shear failures lead to the same overall failure mode. 

Fiber failure pertains only to Option 3 
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CHRISTENSEN STRAIN-BASED FAILURE CRITERION 


New strain-based failure laws are being explored. Criteria used to decide on the failure law to be 
implemented were based on considerations for modeling progressive failure of modem composites under 
impact loading. With the application of thicker section laminate constructions and more ductile fiber- 
matrix systems, it is felt that strain-based theories which clearly delineate between fiber failure and fiber- 
matrix interaction failure are essential for accurate modeling of progressive lamina failures. In addition, 
the failure law should be capable of including three-dimensional strain and stress states typical of thick 
section laminate constructions. 


• New strain-based failure criterion: 

- applicable to more ductile lamina 

- based on Christensen’s quasi-three-dimensional laminate theory 

- coupled theory, interactions of strain (stress) components 

- clearly delineates between fiber failure and fiber/matrix 
interaction failure 
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CHRISTENSEN STRAIN-BASED FAILURE CRITERION (Continued) 


The failure criterion shown below is a modification of one first proposed by Christensen. In this 
law, ejj is the dilatation (first invariant strain) and e kk is the deviatoric strain. We have added the term 

P ej^, the quadratic deviatoric component which was not included in Christensen's work. This allows us 
to have different transverse failure strains in tension and compression, a necessary feature to reproduce 
failure in most lamina constructions. 

• CHRISTENSEN FAILURE CRITERION: 

- Fiber failure 


• Fiber/Matrix Interaction Failure 

“ e m + P4c + e H < 5r 

a, P, Kare determined by simultaneously fitting to three failure states 
usually shear in the 1-2 plane, and positive and negative stress transverse 
to the fiber direction. 
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INITIAL STUDIES OF PROGRESSIVE FAILURE OF 
COMPOSITE AIRCRAFT COMPONENTS 


The figure below summarizes the initial capabilities developed to simulate the impact behavior of 
composite aircraft structural component. 


• Classical Composite Laminate Theory 

- composite laminate can be built up layer-by-layer by specifying 
materials moduli, angle of orientation and stacking sequence of 
each ply 

- each layer can be linear elastic or elasto-plastic 

- maximum strain failure criteria for continuous filament reinforced 
composite material 

• Implemented into a three-node DKT triangular element in DYCAST code 

• Used to simulate NASA drop test of Z-section graphite epoxy frame 



NASA DROP TEST OF A Z-SECTION GRAPHITE-EPOXY 

CIRCULAR FRAME 


This test program is described in "Drop Test and Analysis of Six-Foot Diameter Graphite-Epoxy 
Frames," by Richard L. Boitnott and Huey D. Carden, presented at the AHS National Specialists' Meeting 
on Crashworthy Design of Rotorcraft, April 7-9, 1986. The key features of the test are summarized 
below. The Z-section frame was constructed from four 90-degree sections that were attached by splice 
plates. The laminated composite material was constructed of graphite-epoxy fabric. Another feature 
which added to the computational difficulties was a rear steel restraining plate and forward plexiglass plate. 
The purpose of the plates was to restrain lateral motion. Consequently, contact conditions had to be 
described between these plates and the forward and rear face of the Z-section flanges. 


• Five graphite-epoxy frames were dropped onto a concrete floor 

• Tests were performed at NASA LRC Impact and Dynamics Branch under 
the supervision of Huey Carden and Richard Boitnott 

• Complete frames were fabricated from four 90-degree sections and 
joined with splice plates 

• Z-section cross-sections typical of fuselage structure of advanced 
composite transport aircraft 

• Lateral motions were restrained by a plexiglass and a steel plate 


• Impact was at a splice plate 

• Impact speed was 27.5 FPS 

• Twenty pounds of added weight at the floor beam 
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This figure shows the general configuration of the frame. The frame was modeled with DKT 
triangular finite elements with the material characteristics described above. A number of models were 
investigated. The final model used eight elements through the web at the point of contact. The model 
progressively was made coarser away from the contact point. Beam elements were used in the upper 
quadrant. The final model of one-half of the frame had 867 nodes, 2108 elements, and 5133 degrees of 
freedom. 


Splice Plate 



Gage 


/ 



Frame 

Cross Section 
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PROGRESSIVE FAILURE OF LAMINATED COMPOSITES 


This figure shows a strain trace at the impact point at the inner flange. The location of the strain 
gage is shown on the previous figure. The comparison indicates that the essential features of the response 
were predicted. The time and location and progression of failure also approximated the simulation of a 
composite aircraft structure. However, it should be pointed out that the test involved a number of features 
that increased the difficulty of the analysis so that the primary behavior was not composite failure alone. 
These are: contact between the frame and the face plates, unknown friction coefficient, stiffness of the 
splice plates. Consequently, further analysis and tests on simpler components must be performed in order 
to investigate the use of composite failure criteria in a crash simulation. 


Graphite-Epoxy Fuselage Frame Drop Test 


0.008 
0.006 
0.004 

STRAIN 

0.002 

0.000 
- 0.002 
- 0.004 

0 2 4 6 8 

TIME, ms 



Splice Plate 




] 


Frame 

Cross Section 


86 




STRUCTURAL IMPACT OF A/C COMPOSITES 


This figure outlines areas of continued research. 


* Improved composite failure criteria 

- New, fully interactive strain-based Christensen lamina 
failure criterion to be tested in DYCAST 

- Distinction between fiber and matrix failure modes maintained 

- Evaluate by simulating NASA experiments 




NONLINEAR COMPOSITE BEAM ELEMENT 


A curved open-section composite beam element has been developed. Each flange and web can be 
described as a distinct laminate construction. 


* Nonlinear, composite finite element developed for DYCAST Code 

- Applicable to curved, thin-walled beam/columns with open 
cross-section typical of aircraft frames and stiffeners 

• Each individual flange and web can be a distinct laminate 
composite construction 
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NONLINEAR COMPOSITE BEAM ELEMENT 


T tu; . F , uture work , wi1 .! involve applying failure criterion to each ply of the laminate web or flanges. 

Ultimately, a more detailed global/local representation must be developed to account for local three- 
dimensional effects at the point of impact. 


• Material behavior 
- Laminate composite 


• Failure criteria can be applied to each ply of flange/web defining 
progressive failure during crash loading 
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STUDY OF THE EFFECT OF SHEAR STRESS/STRAIN 
NONLINEARITY ON IMPACT BEHAVIOR 

The next four figures illustrate the influence of material nonlinear behavior and failure mode on the 
dynamic buckling and collapse behavior of composite shells. To examine t h,s > we consider 3 
shell with diaphram ends constructed of typical graphite epoxy lamina with a [0,90,45 r 45]s construction 
The shell is subjected to a hydrostatic step wave pressure loading. Under crash conditions, one wo .^ d 
to have the energy absorption mechanisms for the structure contained in localized failure ^eswith the 
remainder of the structure experiencing little catastrophic damage. To examine the role of nonlinear 
material behavior, we compared the effects of composite failure and shear stress-strain nonlinearity. The 
composite failure law used is the modified Christensen law described previously and material nonlinearity 
was treated by the simple cubic shear stress-strain law due to Tsai and Hahn shown below. 


• Nonlinear shear stress strain behavior can significantly alter 
the dynamic buckling behavior of composites structures 

• Example: 

- cylindrical shell with diaphram ends 

- subjected to a step hydrostatic pressure loading 

• Used Tsai-Hahn stress strain relation 

Yl2 = [^12 1 + ^66 X 12 1 T 12 

G l2 = 0.71 x 10 6 
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STUDY OF THE EFFECT OF SHEAR STRESS/STRAIN NONLINEARITY 
ON IMPACT BEHAVIOR (Continued) 


The effect of shear stress nonlinearity is summarized below and shown on the next visual. 


• Under hydrostatic, step wave loading, the dynamic collapse occurs 
in a mode dominated by the first bifurcation buckling mode of the 
linearly elastic structure if the degree of shear nonlinearity is small 


* For large shear nonlinearities, the collapse resembles a plastic 
instability mode rather than a bifurcation buckling mode 


• Energy absorption and force transmittal characteristics for these 
two modes are quite different and may lead to different crash 
behavior 
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STUDY OF THE EFFECT OF SHEAR STRESS/STRAIN NONLINEARITY 
ON IMPACT BEHAVIOR (Continued) 


Two sets of results are shown below. On the left are results for essentially linearly elastic response 
with the parameter S 6 6 in the Tsai-Hahn relation set equal to 0.5 x 10 13 . The set on the right increases the 

level of nonlinearity by setting S66 to 0.5 x 10 10 . As can be observed, a transition occurs with the 
increase in shear nonlinearity from a traditional buckling mode as shown on the left to a mode that is more 
characteristic of a general collapse as shown on the right. The general collapse mode in which there is a 
larger observed region of high "failure" strains around the circumference of the shell is far more damaging 
to the shell. These initial considerations indicate that material characterization and an understanding of the 
material behavior during impact conditions is an important consideration in predicting structural response 
to a crash loading situation. 



S 66 = 0.5 x 10 13 


S w = 0.5 x 10 10 
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EFFECT OF FAILURE STRAIN ON MODE OF DEFORMATION 


This figure shows the effect that the mode of failure has on the structural response. The figure on 
the right shows the response of the shell when published values of fiber failure are used in the analysis. 
Once again the modified Christensen failure law is used. In this case, initial matrix failure occurs which 
quickly spreads through the cross-section near the center of the shell. The observed failure (shown on the 
right) is similar to a collapse mode shown previously. If the fiber failure strain is reduced so that fiber 
failure occurs first, then the failure mode is more predominately like a buckling mode as seen in the figure 
on the left. 

In the designing of aircraft structures for impact, it is crucial that the dynamic collapse behavior of 
the primary structural components be in modes not highly prone to distributing material failure. As has 
been shown, this requires a thorough understanding of the material behavior when subjected to impact 
loadings. 



FIBER FAILURE 


MATRIX FAILURE 


SUMMARY 


There are four points made below. Before we can simulate the behavior of composite structures 
due to impact loading in order to evaluate crashworthiness, the material and laminate phenomena must be 
identified and understood. Carefully controlled test programs must be designed in order to achieve these 
goals. Because of the complexity of phenomena and the current status of computational techniques, there 
still is a "trade-off between detailed finite element modeling and "hybrid" modeling of energy absorbing 
structures constructed of composite materials. This may also be true for metallic aircraft structures. 


• It is essential to characterize the progressive failure behavior 
of laminated composite materials subjected to impact loadings 
and to implement this capability into finite element simulation 
codes 


• Carefully controlled tests of composite components must be 
performed to better understand modes of failure, energy 
absorbing capabilities and to provide data for correlation with 
crash simulation codes 


• There still may be a "trade-off" between detailed finite element 
and hybrid modeling of laminated composite material in some 
sections of the structure 


• It is essential that crashworthiness be incorporated into the 
structural design for composite materials 
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INTRODUCTION 


During the 1960's over 30 full-scale aircraft crash tests were conducted by the Flight Safety 
Foundation under contract to the Aviation Applied Technology Directorate (AATD) of the U.S. Army 
Aviation Systems Command (AVSCOM). The purpose of these tests were to conduct crash injury 
investigations that would provide a basis for the formulation of sound crash resistance design criteria for 
light fixed-wing and rotary wing aircraft! This resulted in the Crash Survival Design Criteria Designer's 
Guide which was first published in 1967 and has been revised numerous times, the last being in 1989 
(Ref. 1). Full-scale aircraft crash testing is an expensive way to investigate structural deformations of 
occupied spaces and to determine the decelerative loadings experienced by occupants in a crash. This gave 
initial impetus to the U.S. Army to develop analytical methods to predict the dynamic response of aircraft 
structures in a crash. It was believed that such analytical tools could be very useful in the preliminary 
design stage of a new helicopter system which is required to demonstrate a level of crash resistance and 
had to be more cost effective than full-scale crash tests or numerous component design support tests. 

From an economic point of view, it is more efficient to optimize for the incorporation of crash resistance 
features early in the design stage. However, during preliminary design it is doubtful if sufficient design 
details, which influence the exact plastic deformation shape of structural elements, will be available. The 
availability of simple procedures to predict energy absorption and load-deformation characteristics will 
allow the designer to initiate valuable cost, weight and geometry tradeoff studies. The development of 
these procedures will require some testing of typical specimens. This testing should, as a minimum, 
verify the validity of proposed procedures for providing pertinent nonlinear load-deformation data. It was 
hoped that through the use of these analytical models, the designer could optimize aircraft design for crash 
resistance from both a weight and cost increment standpoint, thus enhancing the acceptance of the design 
criteria for crash resistance. 
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SYSTEMS APPROACH TO CRASH RESISTANCE 


For maximum effectiveness, design for crash resistance dictates that a total systems approach be 
used and that the designer consider survivability issues in the same light as other key design considerations 
such as weight, load factor, and fatigue life during the initial design phase of the helicopter. Figure 1 
depicts the system’s approach required relative to management of the crash energy for occupant survival 
for the vertical crash design condition. The crash G loads must be brought to within human tolerance 
limits in a controlled manner to prevent injury to the occupants. This can be accomplished by using the 
landing gear, floor structure, and seat to progressively absorb most of the crash energy dunng the crash 
sequence. Thus, the occupant is slowed down in a controlled manner by stroking/failing the landing gear, 
crushing the floor structure, and stroking the seat at a predetermined load before being subjected to the 
crash pulse which by then has been reduced to within human tolerance limits. . In addition, the large mass 
items such as the overhead gearbox are arrested by stroking/failing of the landing gear or fuselage 
structure, and in some cases, by stroking of the gearbox within its mounts. In this example, assuming that 
the landing gear has been designed to meet the minimum requirements of MIL-STD-1290A, i.e., 20 tt/sec, 
the fuselage would be decelerated to approximately 37 ft/sec at the time of contact with the surface. The 
Army’s most recent helicopters, the UH-60 Black Hawk and AH-64 Apache, are both designed generally 
in accordance with the requirements of MIL-STD-1290A. 


• LANDING GEAR 

• FUSELAGE STRUCTURE 


• SEATS 

• OTHER 
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Figure 1 - Energy Management System 
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MODEL KRASH DEVELOPMENT 


Figure 2 depicts the chronology of the development of model KRASH which is the most 
commonly used computer analysis of the dynamic response of aircraft structure during a crash impact. 
Model KRASH represents the structure with beam elements, crushaWe springs and lumped masses. It is 
intended to provide designers with simplified techniques with which to perform crashworthiness studies 
during the aircraft preliminary design phase. 

In any crash resistant aircraft design a total systems approach must be employed to determine the 
most effective mix of energy attenuation from the landing gear, airframe structure and stroking seat. This 
is where program KRASH can allow you to quickly assess the relative effects of different energy 
attenuating component mixes, thus pointing the way to an optimized system design for crash resistance. 


• DYNAMIC SCIENCE/US ARMY 

• LOCKHEED/US ARMY 

• LOCKHEED/FA A 

• DRI/FAA/NAVY 

• US ARMY R&D PROGRAMS 

• IR&D EFFORTS 

Figure 2 - Model KRASH Development 
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1972-1974 
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1990-1991 

1975-1992 

1980-1992 
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KRASH MODEL CORRELATION 


As part of the development of KRASH, full-scale UH-1 cabin floor sections were crash tested to 
determine the load-deflection characteristics of the subfloor springs. In addition, a full-scale UH-1 was 
crash tested to help develop and validate the KRASH model (Pig. 3). Later, each time an aircraft was to 
be crash tested by the U.S. Army, whatever the reason, the aircraft would be modeled for the planned 
impact conditions using KRASH before the test and then the results would be compared to the test data 
and post test structure deformation. Using the actual impact conditions, if different, along with test 
results, the KRASH program would be exercised to fit the actual test and result. This empirical approach 
to improving the model or to just better understand influences of various structural concepts was used in 
full-scale testing of both the Bell and Sikorsky Advanced Composite Aircraft Program (ACAP) 
helicopters. The latter two tests were unique in that they provided the utility of model KRASH for 
composite structures. 


"KRASH" DYNAMIC MODEL 
OF UH-1 AIRFRAME STRUCTURE 



INTERNAL BEAM ELEMENT (Typical) 
EXTERNAL CRUSHABLE SPRING (Typical) 
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MASS (Typical) 
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Figure 3 - KRASH Model Correlation 
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DESIGN SUPPORT TESTING TO SUPPORT ANALYSIS 


KRASH is being used worldwide and has proven to be a useful tool in the preliminary design 
process for crash resistance in new rotorcraft. However, in preliminary design, energy absorbing 
subfloor spring constants are assumed based upon elemental structural test data, sometimes not much more 
than coupon specimen data. All too often, once the preliminary design advances into detail design and 
some elemental design support testing is conducted, it is found that the structural concept does not perform 
correctly and that the spring constants used initially are very difficult if not impossible to obtain within the 
structure weight allocated. This leads the designer into a trial and error design support test effort (see Fig. 
4) with associated high costs, to obtain the elusive good crash initiators, good energy attenuation and a 
nice rectangular load-deflection relationship. Since this seems to be the case almost all the time, especially 
with composite structures, it seems that if we could develop a database from all this trial and error testing 
and make it available to all designers, we could significantly reduce the duplication of test effort that is 
currently going on with all kinds of test specimens. Furthermore, the analyst could use these test data as 
an empirical base to develop finite element models of the structural elements, thus reducing design support 
testing, though ultimately the final design would still have to be tested under realistic impact conditions. 
The sharing of data and increased use of analytical models will permit crash resistance component design 
formulation in less time at far less cost. 



Figure 4 - Analysis/Test Logic Path 
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SOM-LA 


Program SOM-LA (Seat/Occupant Model-Light Aircraft) has been developed for use in evaluating 
the crashworthiness of aircraft seats and restraint systems. It combines a three-dimensional dynamic 
model of the human body with a finite element model of the seat structure. It is intended to provide the 
design engineer a tool with which he can analyze the structural elements of the seat as well as evaluate the 
dynamic response of the occupant during a crash. 

The occupant model consists of 12 masses that represent the upper and lower torso, neck, head, 
and two segments for each of the arms and legs. An optional model of die human body includes beam 
elements in the spine and neck, but is restricted to two-dimensional motion. 

External forces are applied to the occupant by the cushions, floor and restraint system. Interface 
between the occupant and seat is provided by the seat bottom cushion, back cushion, and an optional 
headrest. The restraint system can consist of a lap belt alone or combined with a single shoulder belt, over 
either shoulder, or a double-strap shoulder harness. A lap belt dedown strap, or negative G strap, can also 
be included. Each component of the restraint system can be attached to either the seat or the aircraft 
structure. SOM-LA is summarized in Fig. 5. 


• FINITE ELEMENT MODEL OF SEAT 

• DYNAMIC MODEL OF SO™ PERCENTILE HUMAN MODEL AND 
ANTHROPOMORPHIC DUMMY 

• GIVES STRUCTURAL BEHAVIOR OF SEATING AND RESTRAINT SYSTEMS 
UNDER TRANSIENT DYNAMIC LOADING CONDITIONS 

• CAPABLE OF PREDICTING SEAT STROKE, OCCUPANT MOTIONS, 
OCCUPANT-FLOOR, OCCUPANT-CUSHION, AND OCCUPANT-RESTRAINT 
FORCES 

• CAN ACCEPT ANTHROPOMETRY INPUT 

• CANNOT ACCEPT FORCES OF OCCUPANT CONTACT WITH SURROUNDING 
STRUCTURE 

Figure 5 - Seat Occupant Model/Light Aircraft (SOM-LA) 
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ARTICULATED TOTAL BODY MODEL 


The Articulated Total Body (ATB) Model is primarily designed to evaluate the three-dimensional 
dynamic response of a system of rigid bodies when subjected to a dynamic environment consisting of 
applied forces and interactive contact forces. Although the ATB Model was originally developed to model 
the dynamic response of crash dummies and, with later modifications, the response of the human, the 
ATB Model is quite general in nature and can be used to simulate a wide range of physical problems that 
can be approximated as a system of connected or free rigid bodies. 

'Hie approach used in the ATB Model to model the human or manikin body (the "body" in the ATB 
Model simulation) is to consider the body as being segmented into individual rigid bodies (the "segments" 
in the ATB Model) each having the mass of the body between body joints or, in the case of single-jointed 
segments, such as the foot, distal to the joint. An example would be the left upper arm segment, which 
represents the mass of the body between the shoulder joint and elbow joint. Segments are assigned mass 
and moments-of-inertia and joined at locations representing the physical joints of the human body, such as 
the shoulder joint or the knee joint. For the ATB Model the Generator of Body Data (GEBOD) is a source 
of anthropomorphic data for the zero to 100th percentile male, female, infant, child and dummy. These 
data include body masses, moment-of-inertia, and c.g.'s. 

A personal computer version of ATB is named DYNAMAN which along with the ATB model has 
been useful in R&D programs to delethalize the helicopter cockpit. The ATB model is summarized in Fig. 


• PREDICTS GROSS HUMAN BODY RESPONSE TO VARIOUS 
DYNAMIC ENVIRONMENTS 

• CAN ACCEPT INPUT OF STROKING SEATS AND RESTRAINTS 

• GEOMETRIC BODY MODELER (GEBOD) IS A SOURCE OF A 
WIDE RANGE OF ANTHROPOMORPHIC DATA INPUT 

• CAN PREDICT CONTACT FORCES ON OCCUPANT 

Figure 6 - Articulated Total Body Model 
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INTRODUCTION 


Modem military aircraft transparency systems, windshields and canopies, are complex 
systems which must meet a large and rapidly growing number of requirements. As illustrated in 
Fig. 1, many of these transparency system requirements are conflicting-^presenting difficult 
balances which must be achieved. One example of a challenging requirements balance or trade is 
shaping for stealth versus aircrew vision. 

The large number of requirements involved may be grouped in a variety of areas including 
man-machine interface; structural integration with the airframe; combat hazards; environmental 
exposures; and supportability. Some individual requirements by themselves pose very difficult, 
severely nonlinear analysis problems. One such complex problem is that associated with the 
dynamic structural response resulting from high energy bird impact. 

OBJECTIVE: 


MISSION SUPPORTABILITY 

CRITICAL AS A DELIVERABLE 

PERFORMANCE FEATURE 



Figure 1 - Balance between performance and supportability 
requirements for an aircraft transparency system. 


1 Brockman, R. A.; and Held, T. W. : Explicit Finite Element Method 
for Transparency Impact Analysis r WL-TR-91-3006, 1991. 
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NONLINEAR FINITE ELEMENT ANALYSIS OF AIRCRAFT 
TRANSPARENCY BIRD IMPACT 


Impact phenomena encompass a broad range of structural behavior and response times, 
which depend upon the stiffness, strength, mass, geometry, velocities, and failure characteristics 
of the bodies involved. Soft body impacts, such as transparency bird impacts, are unusual: while 
the response is often highly nonlinear, critical features of the response may occur either at early 
times or long (milliseconds) after the impact is finished as illustrated in Fig. 2. For over ten years, 
implicit solution techniques with isoparametric solid finite element technology (Ref. 1) have been 
used successfully to analyze aircraft transparency bird impact response (Refs. 2-5). An impact 
solution may be dominated by complicated contact conditions which preclude the use of large time 
steps, so that the advantages of an implicit solution are lost. Practical transparency analysis remains 
a time-consuming and laborious process, and in some circumstances the present inventory of 
analysis tools may not be optimal. 



Figure 2 - Dynamic bird impact response of F-16 fighter aircraft 
prototype canopy design at 20 msec. 
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EXPLICIT FINITE ELEMENT METHOD FOR AIRCRAFT 
TRANSPARENCY BIRD IMPACT 


This paper outlines the development of new computational techniques for analyzing 
structural response to high-speed impact. The key improvements in the new technique are listed in 
Fig. 3. The analytical technique discussed is an explicit finite element method of the type used 
widely for the numerical solution of shock and wave propagation problems. The explicit family of 
time integration algorithms is attractive because it is readily adapted to high performance on the 
current generation of supercomputers, which combine parallel or pipeline processors, moderate 
amounts of high-speed memory, and relatively slow disk performance. An added benefit is the 
ability to implement more detailed material and failure models. The particular implementation 
discussed here is a computer code called X3D. X3D is an explicit, three-dimensional finite element 
program intended for use in solving impact, wave propagation, and other short-duration problems 
in structural dynamics. 


• Soft-body impact loads: the bird appears explicitly in the finite element 
model, so that ad hoc estimates of the impact loading distribution are 
unnecessary 


• Material modeling: the material models include strain rate sensitivity and 
failure 


• Layered shells: multilayered constructions, including those with soft 
interlayers, can be modeled using a single layer of surface elements 


Figure 3 - Key improvements offered by explicit finite element methods 
for nonlinear dynamic aircraft transparency bird impact. 
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X3D EXPLICIT THREE-DIMENSIONAL FINITE ELEMENT PROGRAM 
FOR SHORT-DURATION STRUCTURAL DYNAMIC PROBLEMS 


X3D contains two types of finite elements: solids and plates. The solid elements are an 
eight-node hexahedron, based on a mean-stress approximation with anti-hourglass stabilization 
(Ref. 6); and a four-node tetrahedron. The eight-node solid hexahedron element is illustrated in 
Fig. 4. The solid hexahedral finite element uses a displacement and velocity approximation based 
upon trilinear polynomials; that is, the element's displacement and velocity components each vary 
linearly along each edge of the element. In addition, the stress components are computed from a 
mean stress approximation using only the mean velocity gradient for the element (Ref. 6). This 
measure is desirable to maintain good element performance, and also reduces the effort required for 
element computations. However, the resulting element is a constant stress element, and therefore a 
generous number may be necessary for accurate modeling. In particular, a single layer of these 
solids is incapable of developing a bending moment. The material model used for solids consists of 
a polynomial equation of state coupled with a von Mises plasticity model, a simple power-law 
correction for strain rate sensitivity, and a failure criterion based upon the ultimate stress. 



W 


Nodal degrees 
of freedom 



U 


Figure 4 - Eight-node solid hexahedron X3D element. 
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X3D SOLID HEXAHEDRAL AND TETRAHEDRAL ELEMENTS 


■ . because of the mean stress approximation, certain modes of deformation exist for the 

hexahedron which are stress-free but do not represent rigid body motions. These hourglassing 
deformation modes correspond to linearly varying stress patterns which are not detected by the 
mean stress approximation as shown in Fig. 5. To stabilize these potentially unstable motions an 

° T 1 ? 03 ^ -° reSiSt the hourglass motions trough internal damping 
i 7 ' The te ? ah f lral s 2 ll , d 1S a constant-strain, constant-stress element based upon fully 
linear displacement and velocity field approximations. The element is quite similar to the ? y 
hexahedron, but does not use an anti-hourglass viscosity. The twelve degrees of freedom for the 

nnSili H 3 f ture l ? e S1X ngld ‘ body I I lotions die six uniform strain/stress modes, so that no 
... . e 0r ^ iatl0 . n P atterns exist f° r individual elements. The tetrahedron is included in X3D for 
m f )ft ‘ b ? dy ln ? pact rnodeling. Since the element has no unstable modes, it can be used to 
u 1S H 0rtl0nS wltho 1 ut causing numerical problems. The tetrahedron is used to 
model the bird in bird impact simulations, using an equation of state typical of water a very low 
strength deviatonc model and an ultimate failure strain of about 5 (500%). The tetrahedron is 
implemented as a five-node element, the fifth node coinciding with the first. This artifice serves to 
istinguish the four-node tetrahedron from the four-node quadrilateral plate element during input. 


HOURGLASS DETORMATION PATTERNS TOR 8-NOOE SOUL) 
DISPLACEMENT DIRECTION -* ► 





Figure 5 - Hourglass deformation patterns for solid element. 
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MATERIAL MODELS 


The material constitutive relationships used for both the solid and plate finite elements consist of 
. Son and a bulk(pressure-volume) model The stress tensor is composed of 

„ hvdmstatic or pressure stress, and a deviatoric stress tensor which is independent o 
pressure. These two contributions to the stress tensor are determined independently in the material 
model by a deviatoric stress model and a mechanical equation of state. The deviatoric material 
model uLd for solids is a rate-dependent, isotropic, hypoelastic theory appropriate for moderate to 
IZe def^do^s The parameter which define the material’s deviatoric! behavior are shown in 
Fie 8 An experimental feature provides an isotropic Newtonian fluid model for the three-^ 
dimensional Sid elements for potential use in hydrodynamic impact modeling. The bulk behavior 
is described in polynomial form as for a solid, while the deviatoric stress is related linearly to the 
rate of deformation In the plate element, the elastic-plastic material model is slightly mo 
SmplicaShan for three-dimensional solids because of the zero normal stress _ constraint ^During 
the plasticity calculation, it is necessary to determine a final state of stress which notonly lies on 
S e SfsurfS but which satisfies the condition for the nonnal stress to be zero. The ttevuuonc 
iriodel and the bulk model (equation of state) are not entirely independent, and must be solved 
simultaneously with the normal stress constraint. 


• Linear shear modulus 

• Quasi-static yield stress 

• Rate sensitivity scale factor 

• Rate sensitivity exponent 

• Hardening modulus 

• Ultimate stress 


Figure 8 - Parameters for material deviatoric behavior. 
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X3D SOLID HEXAHEDRAL AND TETRAHEDRAL ELEMENTS 


Because of the mean stress approximation, certain modes of deformation exist for the 
hexahedron which are stress-free but do not represent rigid body motions. These hourglassing 
deformation modes correspond to linearly varying stress patterns which are not detected by the 
mean stress approximation as shown in Fig. 5. To stabilize these potentially unstable motions, an 
anti-hourglass viscosity is employed to resist the hourglass motions through internal damping 
forces (Ref. 7). The tetrahedral solid is a constant-strain, constant-stress element based upon fully 
linear displacement and velocity field approximations. The element is quite similar to the 
hexahedron, but does not use an anti-hourglass viscosity. The twelve degrees of freedom for the 
element capture the six rigid-body motions and the six uniform strain/stress modes, so that no 
unstable deformation patterns exist for individual elements. The tetrahedron is included in X3D for 
its utility in soft-body impact modeling. Since the element has no unstable modes, it can be used to 
follow very large distortions without causing numerical problems. The tetrahedron is used to 
model the bird in bird impact simulations, using an equation of state typical of water, a very low 
strength deviatoric model, and an ultimate failure strain of about 5 (500%). The tetrahedron is 
implemented as a five-node element, the fifth node coinciding with the first. This artifice serves to 
distinguish the four-node tetrahedron from the four-node quadrilateral plate element during input. 


HOURGLASS DEFORMATION PATTERNS TOR B NODE SOLID 
DISPLACEMENT DIRECTION * 





Figure 5 - Hourglass deformation patterns for solid element 




X3D PLATE AND SHELL ELEMENT 


The plate and shell element in X3D is a four-node quadrilateral based upon a Mindlin- 
Reissner type thick-plate theory. A corotational axis system, which rotates with the element but 
does not deform, is used to simplify the element kinematics. The plate and shell element uses a 
reduced (one-point) Gaussian quadrature, in conjunction with anti-hourglass stabilization 
techniques. An approximate model for layered media is implemented for the element, so that plates 
and shells having layers with large differences in stiffness can be represented effectively using a 
single element in the thickness direction. For each layer of the X3D plate and shell element, the 
material is elastic, perfectly plastic, and obeys plane stress assumptions. Transverse shear stresses 
in the element are uncoupled from the tangential stresses, and follow an elastic constitutive relation. 
The plate and shell element has six degrees of freedom per node as shown in Fig. 6. The 
displacement and rotation components each are interpolated separately, using bilinear polynomials. 
The resulting element is quite similar to that described by Belytschko (Ref. 8). 



Figure 6 - Four-node quadrilateral plate element. 
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X3D PLATE AND SHELL ELEMENT 


Unlike the solid elements, the plate element must be formulated in a local axis system 
because of the differing treatment of the plate thickness from that of the planform directions. A 
corotational coordinate system which rotates with the element is employed, and therefore is 
constructed anew at each time step of the solution based upon the current element geometry. The 
plate element shape functions are formulated entirely in local coordinates. The element calculations 
are performed with respect to the "mean plane" of the element, and corrected as necessary to 
account for out of plane warping of the reference surface. The plate uses a mean-stress 
approximation for its inplane directions, similar to the solid hexahedron. At any thickness station, 
the velocity gradient is evaluated at the centroid of the element, and assumed to be constant 
throughout the element (except through the thickness). To resist unstable motions resulting from 
the assumption of a uniform velocity gradient, the plate element uses a stiffness hourglass control 
scheme (Ref. 6). Other aspects of element design are listed in Fig. 7. 


• Simpson's Rule integration through the thickness 


• Each layer may be a different material and even use a different 
material model 


• Layered constructions with dramatic stiffness characteristics variation 
from layer to layer require special treatment 


• Formulation of lumped mass coefficients relieves stringent time step 
restriction without upsetting convergence (Ref. 9) 


• No inertia is assigned to the "drilling" rotation in the local coordinate 
system 


Figure 7 - X3D plate and shell element features. 



MATERIAL MODELS 


The material constitutive relationships used for both the solid and plate finite elements consist of 
a deviatoric (shear) relation and a bulk (pressure-volume) model. "Die stress tensor is composed of 
a hydrostatic, or pressure, stress, and a deviatoric stress tensor which is independent of the 
pressure. These two contributions to the stress tensor are determined independently in the material 
model by a deviatoric stress model and a mechanical equation of state. The deviatoric material 
model used for solids is a rate-dependent, isotropic, hypoelastic theory appropriate for moderate to 
large deformations. The parameters which define the material's deviatoric behavior are shown in 
Fig. 8. An experimental feature provides an isotropic Newtonian fluid model for the three- 
dimensional solid elements for potential use in hydrodynamic impact modeling. The bulk behavior 
is described in polynomial form as for a solid, while the deviatoric stress is related linearly to the 
rate of deformation. In the plate element, the elastic-plastic material model is slightly more 
complicated than for three-dimensional solids because of the zero normal stress constraint. During 
the plasticity calculation, it is necessary to determine a final state of stress which not only lies on 
the yield surface, but which satisfies the condition for the normal stress to be zero. The deviatoric 
model and the bulk model (equation of state) are not entirely independent, and must be solved 
simultaneously with the normal stress constraint. 


• Linear shear modulus 

• Quasi-static yield stress 

• Rate sensitivity scale factor 

• Rate sensitivity exponent 

• Hardening modulus 

• Ultimate stress 


Figure 8 - Parameters for material deviatoric behavior. 
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LAYERED PLATE AND SHELL MODEL 


X3D provides 3 method of approximation for plates and shells having large stiffness 
variations from layer to layer, such as those for a laminated aircraft transparency system shown in 
Fig. 9. Layered structures of this type often require detailed and expensive models, since 
conventional plate and shell finite elements do not reproduce the correct transverse shear strain 
distributions through the wall thickness. The X3D method requires only a single layer of elements 
having six engineering degrees of freedom per node, regardless of the number of layers in the 
structure. The approximation uses closed-form elasticity solutions to develop transverse shear 
flexibility corrections, which bring this contribution to the energy into line with that caused by pure 
bending, twisting, and extension. For large displacement problems, the technique is applied in 
corotational coordinates. Changes in stiffness caused by plasticity can be accounted for by 
recomputing the flexibility corrections based upon instantaneous moduli. Applied forces in X3D 
may consist of body forces and surface pressure. Kinematic boundary conditions may include 
prescribed nodal displacements, rigid-wall constraints, and contact between specified surfaces 
within the mesh. Initial conditions may be specified for the translational velocity components for all 
or part of the finite element model. 
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TAYLOR CYLINDER SAMPLE ANALYSIS 


The Taylor cylinder experiment, which is used to estimate the mechanical properties of 
metals at high strain rates, involves the normal impact of a cylinder onto a rigid surface. It is a 
common benchmark problem with a well-known solution. An X3D model was prepared for one 
quarter of the cylinder using 1 350 8-node solid elements. Material constants typical of copper were 
used. Purely isotropic strain hardening was assumed, and no ultimate stress was specified (i.e., 
elements could not fail during the solution). Virtually all of the kinetic energy of the cylinder is 
dissipated through plastic deformation within about 80 microseconds. Figure 10 shows a deformed 
mesh plot of the cylinder in its final state. 



|Y 


X Z 

TAYLOR CYLINDER R * 3 2 MM . L » 324 MM DOUBLE SYMMETRY 

DISPLACEMENTS FOR TIME STEP - 8886 / TIME - 8 OOOOOOOE 05 

11/13/90 2253:12 


Figure 10 - Deformed geometry of Taylor Cylinder. 
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TAYLOR CYLINDER SAMPLE ANALYSIS 


Figure 1 1 shows a time history of the cylinder’s length. The analysis was performed in 
8886 time steps, and required 6 hours, 36 minutes on a VAX 8650 computer (about 0.00198 CPU 
seconds per element time step). The same analysis runs in about 40 minutes on a CRAY X-MP 
(.0002 seconds per element time step). Results from the X3D solution compare very well with 
analyses using the DYNA and NIKE codes, as shown below. 


QUANTITY 

X3D 

DYNA2D 

DYNA3D 

NIKE2D 

Final length, mm 

21.47 

21.47 

21.47 

21.47 

Maximum radius, mm 

7.081 

7.127 

7.034 

7.068 

Maximum strain 
at center 

2.95 

3.05 

2.95 

2.97 


CYLINDER LENGTH 



Time, /j, s 

Figure 1 1 - Cylinder length versus time. 
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EXPLOSIVELY LOADED CYLINDRICAL SHELL SAMPLE ANALYSIS 


Marchertas and Belytschko present both computational and experimental results for this problem 
(Ref. 10). A 120 degree cyiindncal panel is loaded by igniting a charge spread over most of the 
surface. In the numerical solution, we represent this impulsive loading by a uniform initial velocity 
along the radius of the shell. A three point integration through the thickness of the shell was used 
with X3D. This is the minimum thickness integration order, and may give a solution which is 
slightly too flexible. Figure 12 shows the geometry of the explosively loaded cylindrical shell. The 
geometric and material parameters for the shell were: 


Radius 2.9375 in. 

Thickness 0.125 in. 

Length 12.56 in. 

Velocity 5,650 in./sec (initial) 


Tensile modulus 
Density 
Yield stress 


10.500.000 psi 
0.0965 Ib/cu.in. 

44.000 psi 



Figure 12 - Geometry of explosively loaded cylindrical shell. 
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EXPLOSIVELY LOADED CYLINDRICAL SHELL SAMPLE ANALYSIS 


Results of the X3D solution, which was performed in 886 time steps, are shown in Fig. 

13. The response mainly involves a flattening of the inner portion of the shell, consisting mostly of 
permanent deformation. The displacements peak at around 0.4 ms, with the largest inward 
displacements approaching half the radius. After this point, there is some elastic recovery (lasting 
about another 0. 1-0.2 ms), but only very small vibration, since most of the energy has been 
dissipated through plastic flow. Displacement histories at selected points agree quite well with 
experimental results. Note that the initial velocity components are directed radially inward, and that 
points on the edges of the loaded region were assigned half the nominal initial velocity to provide 
the correct impulse to the shell. 



Z 

IMPULSIVELY LOA&ED CYLINDER 
DISPLACEMENTS FOR TIME STEP 886 
9-JUL-89 04 49 30 


Figure 13 * Final deformed shape of cylindrical shell. 
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F-16 AIRCRAFT CANOPY BIRD IMPACT SAMPLE ANALYSIS 


The F-16 bubble canopy provides a useful example for validation since the impact response 
involves very large motions, and the coupling between the load distribution and the deformation is 
strong. As a first step in validating the X3D code for bird impact simulation, several analyses of 
centerline impacts were carried out for the original production canopy, a 0.5 in. thick monolithic 
polycarbonate design. This design is capable of defeating 4 lb bird impact at airspeeds up to about 
350 knots. Figure 14 shows the geometry of the transparency and of the projectile, a 4 lb bird 
idealized as a right circular cylinder. The patch outlined around the crown of the canopy and the 
entire bird are covered with contact elements. The canopy model consists of 928 quadrilateral plate 
elements. The bird is represented by 960 tetrahedral solids with equation-of- state coefficients 
typical of water, and very small shear stiffness and strength. The low-strength bird model, used in 
about half of the simulations, produces a pressure-volume response similar to water, and a "battle" 
shear behavior. The ultimate and yield stresses coincide, so that element failure occurs at relatively 
small strains. 



Figure 14 - Contact element grid for bird impact problem. 
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F-16 AIRCRAFT CANOPY BIRD IMPACT SAMPLE ANALYSIS 


A high-strength bird model was used as well, which permits roughly 500% plastic 
(deviatoric) strain before the material is declared failed. The question of which bird model is more 
realistic has not been resolved because so many details are unknown with regard to material 
properties, precise support conditions, and center of impact location. Note that when elements of 
the bird model fail due to large shear distortion, their mass is retained in the problem, and the 
corresponding nodes continue to be used in contact calculations. Therefore, portions of the 
impacting body which have "failed" continue to transfer momentum to the target, but do not 
contribute to the summation of internal forces. In the deformed plot shown in Fig. 15, nodes 
attached to failed elements in the bird model are shown as small circles representing the center of 
mass positions. For the cases considered, the center of impact is at fuselage station 1 12 (measured 
in inches), which is about two feet aft of the forward edge of the canopy. The initial velocity of the 
bird is horizontal and equal to 350 knots (7,094) in./sec) at all nodes. The solution illustrated in 
Fig. 15 employs the low-strength bird. The displacement results are similar to experimentally 
observed values, although the computed deformed shape exhibits larger displacements in the 
forward region. 



F-16 BIRD IMPACT AT FS 112 
DISPLACEMENTS FOR TIME STEP L 
05-JUN-61 04;27:38 


SYMMETRIC, 350 KNOTS 

1731 t TIME * 6 0 0 0 0 0 0 1 E - 03 


M60E ' 0 5 INCH • POLYCARBONATE SV=7140 / SU=16000 


Figure 15 - Deformed geometry of F-16 canopy for low-strength bird 350 knot impact. 
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SUMMARY AND CONCLUSIONS 


An improved analytical capability for soft-body impact simulation has been developed. 
Advances have been made in modeling impact loads, nonlinear materials, and layered wall 
constructions. The explicit approach adopted exploits the strengths of the current generation of 
supercomputer hardware, so that analysis cost and turnaround times are reduced significantly over 
the previous generation of bird impact analysis software. Experience in performing applications 
indicates that the solution methodology is very reliable, requiring minimal user intervention to 
avoid or correct problems with the solution. Implicit methods used in earlier work on these 
problems demand a great deal of user attention for stable, accurate, and convergent results, while 
the explicit technique is relatively trouble-free. The work reported here is a significant step toward 
a reliable capability for design screening and parametric investigation. Figure 16 lists the two 
primary research needs required to complete such a capability. With a modest effort in these areas 
of research need, the techniques and software described can become a truly useful and reliable tool 
for design and evaluation of a new generation of bird-impact resistant aircraft transparency 
systems. 


• Model Validation. Additional comparisons of analytical predictions 
with full-scale impact test data are needed to develop confidence in 
the accuracy of the analysis and knowledge of its limitations. 


• Materials Characterization. The transparency materials in wide use are 
high-polymer compounds with very complex characteristics. Much more 
experimental and analytical work is needed to understand these materials 
adequately and model their behavior faithfully. 


Figure 16 - Research needs for aircraft transparency bird impact, explicit finite 

element analysis methods. 
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VEHICLE CRASHWORTHINESS DESIGN 


Objectives 

To design the vehicle structure for optimum impact energy absorption, 
and to design the restraint system (seatbelts, airbags, bolsters, etc.) for 
optimum occupant protection. 


Approach 

A major part of the impact energy is to be absorbed by the vehicle structure. 
The restraint components will provide protection against the remaining 
crash energy. 

Certain vehicle components are designed to deform under specific types 
and speeds of impact in a desired mode for sound energy management. 

Structural components such as front side rails, rear rails, door 
structure and pillars undergo large amounts of deformation. 

With properly designed geometry and material these components assist 
in mitigating the effects of impact. 
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Frontal Crash FEA Model 

The original and deformed geometry of a frontal barrier impact 
model is shown below. Typically, model size is approximately 25000 
elastic plastic shell and solid elements. Computer runtime for 
100 milliseconds of crush is about 20 hours on CRAY YMP. The analysis 
is performed with a commercially available code, RADIOSS, from 
Mecalog, France. 
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VEHICLE CRASHWORTHINESS DESIGN 


Objectives 

To design the vehicle structure for optimum impact energy absorption, 
and to design the restraint system (seatbelts, airbags, bolsters, etc.) for 
optimum occupant protection. 


Approach 

A major part of the impact energy is to be absorbed by the vehicle structure. 
The restraint components will provide protection against the remaining 
crash energy. 

Certain vehicle components are designed to deform under specific types 
and speeds of impact in a desired mode for sound energy management. 

Structural components such as front side rails, rear rails, door 
structure and pillars undergo large amounts of deformation. 

With properly designed geometry and material these components assist 
in mitigating the effects of impact. 
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Role of CAE 


CAE is playing a vital role in crashworthiness design of automotives. 

Instead of being dependent on numerous prototypes, the design process 
is rapidly becoming CAE guided. 

This new approach allows engineers to examine many alternatives 
in a shorter time period. 

Also, CAE models are complimenting prototype tests. 

Models are used to predict crash performance of: 

Components 
Subassembly, and 
Full-vehicle systems. 


CAE Tools for Crashworthiness Analysis 

Four types of CAE models are generally used in the industry. These are: 

• Concept models (lumped masses and springs) 

• Occupant simulation models (rigid body type) 

• Hybrid models (concept and FEA combined with test data) 

• Detailed finite element models. 

Use of these models depends on the stage of vehicle development and 
accuracy of results desired. 

CPU usage and modeling turnaround time is also a key factor. 
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Finite Element Models 


G mrgntly at Ford Motor Company nonlinear finite element models are 
used to simulate : 

Front Crash 

• Frontal perpendicular barrier impact 

• Frontal angular barrier impact 

• Frontal pole impact 

Side Impact 

• Side impact simulation 
Rear Impact 

• Vehicle rear-end impact with moving rigid barrier 

• Car-to-car rear offset impact 

Roof Crush 

• Quasi-static roof crush of vehicle structure 


Frontal Crash Models 


Applications: 

• Front-end structural design of the vehicle 

• Provide design directions for airbag sensor development 

• Provide input for occupant simulation models 

FEA Model Output: 

• Total vehicle and component collapse pattern 

• Engine compartment energy absorption 

• Front bumper and side rail load-displacement 

• Engine block movement and dash wall intrusion 

• Deceleration pulse in passenger compartment 

• Displacement and velocity histories at critical locations 



Frontal Crash FEA Model 

The original and deformed geometry of a frontal barrier impact 
model is shown below. Typically, model size is approximately 25000 
elastic plastic shell and solid elements. Computer runtime for 
100 milliseconds of crush is about 20 hours on CRAY YMP. The analysis 
is performed with a commercially available code, RADIOSS, from 
Mecalog, France. 
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Frontal Pole Impact Model 


Finite element model of a vehicle impacting a rigid pole at 31 mph is 
shown below. The model has approximately 28000 elements. The 
deformed geometry shows that the pole is wrapped around by the front 
bumper beam and the engine block is pushed rearward. 



ACCELERATION (G'S) DISPLACEMENT (INCHES) 


Comparison of B-Pillar/Rocker Joint Longitudinal 
Displacement with Test Data 



Accelerations at Right Side Radiator 
Support From Model and Test 






Comparison of Predicted Velocity at Left 
Side Radiator Support with Test Data 



Comparison of Predicted Velocity at Right 
Shock Tower with Test Data 





Front Crash System Model 


Integrated system model consisting of vehicle, occupant, airbag, 
steering column systems, and other interior environment to predict 
both structural and dummy responses from the same run. 
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Vehicle Rear impact Simulation 


The car-to-car offset impact model is developed to predict strength and 
energy absorption of vehicle rear-end (rear rails, rear bumper system, 
lower back, quarter panel and rear floor). The main tasks are to minimize 
fuel tank deformation under severe impact, predict relative motion of 
bullet and target vehicles and displacement, velocity and deceleration 
histories all over the vehicles. 
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FEA Dummy-Structure Side Impact System Model 

The side impact system model consists of a moving deformable honeycomb 
barrier, side impact dummy model and the vehicle structure. Approximate 
model size is 1 6000 elements and CPU needed on CRAY YMP is about 1 5 
hours. 

The detailed FEA model predicts strength and energy absorption by side 
structure (door, pillars, quarter panels) and energy absorption by trim 
and bolster as well as front and rear seat dummy responses (pelvic, 
spine and rib accelerations). 
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Displacement (in) Velocity (mph) 


Vehicle Structure Response From Side impact Model 


Comparisons of model predictions of B-pillar lateral velocity and 
displacement histories with corresponding test data. 




(b) B-pillar Lateral Displacement 



Acceleration (G’s) Acceleration (G’s) 


Occupant Responses From Side Impact System Model 


Comparisons of front dummy pelvis acceleration, T12 acceleration 
and upper and lower rib accelerations from model and tests. 
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INTRODUCTION 


^ The automotive industry has used computational methods for crashworthiness since the early 
1970 s. These methods have ranged from simple lumped parameter models to full finite element models. 
Tne emergence of the full finite element models in the mid 1980's has significantly altered the research 
direction. However, there remains a need for both simple, rapid modeling methods and complex detailed 
methods. This paper will discuss some directions for continuing research. 
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OCCUPANT MODELS 


Although most of the published work has focused on the structural models, some work is starting 
to appear on interior components and the occupant. Typically, the occupant modeled is the mechanical 
surrogate used in tests rather than the actual human. In these cases the accuracy of the finite element model 
is compared to the performance of the mechanical surrogate, called the anthropomorphic test dummy 
(ATD) rather than to the human response. The performance requirements for the ATD are specified based 
on bio-mechanical modeling of the human. Federal mandated standards are based on mechanical quantities 
(acceleration, force) that are measured on the ATD. Thus, in the figure below the human cadaver results 
were used to establish a corridor of behavior in which the Hybrid III ATD chest should perform. 


F(t) 




Figure 3. ATD thoracic response for pendulum impact (Ref. 2). 
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INTRODUCTION 


The automotive industry has used computational methods for crashworthiness since the early 
1970's. These methods have ranged from simple lumped parameter models to full finite element models. 
The emergence of the full finite element models in the mid 1980's has significantly altered the research 
direction. However, there remains a need for both simple, rapid modeling methods and complex detailed 
methods. This paper will discuss some directions for continuing research. 
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COMPUTATIONAL METHODS FOR CRASHWORTHINESS 


Although the design problem in automotive safety is coupled, traditionally automotive 
crashworthiness analysis has been separated into methods that address the structural performance and 
methods that address the occupant performance. The coupling is usually handled by using the results of 
the structural analysis, typically passenger compartment accelerations and some panel motions, as input 
into the occupant analysis. In addition, a wide range of modeling methods are used ranging from simple 
lumped parameter models to full nonlinear finite element models and including virtually all possible 
intermediate hybridizations. The simple models have the advantages of ease of generation and 
interpretation with some sacrifice of accuracy, whereas the finite element methods have the potential to 
express our best knowledge of mechanics but at the expense of model creation and computational time. 



Lumped parameter occupant 
model with air cushion 




^ a. 

Lumped mass structural model 


Figure 1. Typical computational models. 
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USE OF NONLINEAR FINITE ELEMENT METHODS 


Since tjjf thrust of this workshop is to focus on new directions in computational methods this 

c 9 nc . entrate . 0n fo* fl J} ite element methods. By now virtually every automotive company 
has published simulations using one of the many commercially available nonlinear finite element 

on«w!™ S ' T ^ ese 4! mU atl0 | 1S rea t chin g significant complexity, showing front, side, rear, car-to-car and 
angled impacts. The typical results that are shown are a deformed mesh plot and some overall measure of 
accuracy, typically a passenger compartment deceleration. 


Front impact 




Oblique impact 




Side impact 




Rear impact 
(car-to-car) 



Figure 2. Typical nonlinear finite element results (Ref. 1). 
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OCCUPANT MODELS 


Although most of the published work has focused on the structural models, some work is starting 
to appear on interior components and the occupant. Typically, the occupant modeled is the mechanical 
surrogate used in tests rather than the actual human. In these cases the accuracy of the finite element model 
is compared to the performance of the mechanical surrogate, called the anthropomorphic test dummy 
(ATD) P rather than to the human response. The performance requirements for the ATD are specified based 
on bio-mechanical modeling of the human. Federal mandated standards are based on mechanical quantities 
(acceleration, force) that are measured on the ATD. Thus, in the figure below the human cadaver resu 
were used to establish a corridor of behavior in which the Hybrid III ATD chest should perform. 



Figure 3. ATD thoracic response for pendulum impact (Ref. 2). 
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SO WHAT ARE THE PROBLEMS? 


Reading the literature, one would conclude that nonlinear finite element methods are a well used 
and understood technology and that there is relatively little research work that needs to be pursued. We 
believe that there are four areas that need to be examined. First, the published work tends to mask a real 
difficulty in capturing local details of response that can be significant in crash performance. There is a 
dearth of appropriate failure models of any except the traditional isotropic materials. The issue of 
implementation for these techniques on the next generation of parallel computers is important. Finally, 
how these models are integrated in a rapidly shrinking design cycle must be addressed. 


1. ACCURACY OF LOCAL DETAILS 

2. MODELS FOR NON-TRADITIONAL MATERIALS 

3. MASSIVELY PARALLEL IMPLEMENTATIONS 

4. DESIGN 


Figure 4. Issues to be addressed. 
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FIDELITY OF LOCAL MODELS 


The literature rarely contains the efforts that are expended to produce high quality local results that 
ultimately lead to high quality global results. In this figure two models of an engine cradle are shown. 

The engine cradle is used to attach the engine to the front rails in a front wheel dnve vehicle. The 
difference in these two models is the inclusion of the small lightening holes which are approximately 1 cm. 
diameter. The subsequent difference in geometric failure modes defines how the engine will move in the 
crash and may lead to a significant difference in occupant behavior. 



Figure 5. Two detailed models of an engine cradle (Ref. 3). 
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FIDELITY OF LOCAL MODELS 


A second example is the rail shown below. The difference in the two models is the level of mesh 
density. Clearly the failure modes are different. While it is easy to dismiss these two examples as trivial 
examples of modeling techniques, virtually every project that we are aware of has encountered similar 
difficulties that could not be resolved without careful examination of test results and repeated remodeling. 
To become an effective design tool we need to be able to better understand the complexity of modeling 
required to make design decisions. 



Figure 6. Influence of mesh density on front rail response (Ref. 3). 


MATERIAL MODELS 


As indicated earlier the current suite of effective material models available in the nonlinear finite 
element programs is somewhat limited to isotropic materials. Although there are current efforts to model 
non-traditional materials, our success rate for modeling nonlinear and failure behavior of non-traditional 
materials for design is not strong. In particular, we see a need to model fabric (airbag), composites 
(chopped and continuous), honeycombs, and foams. 


NONTRADITIONAL MATERIALS 


AIRBAGS 

CHOPPED AND CONTINUOUS FIBER COMPOSITES 

HONEYCOMB 

FOAMS 


Figure 7. Non-traditional materials. 
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COMPUTATIONAL EFFICIENCY 


There is a widespread belief that the next big leap in computational efficiency will come from the 
massively parallel computer architecture. There is also a belief that this will require fundamental changes 
to the current group of codes that are being developed. As the figure shows, for example, significant 
speedups can be obtained if contact algorithms are properly implemented. If the changes needed to 
implement these massively parallel versions of existing codes are as great as presently thought, it is 
possibly a project for which a cooperative approach would be appropriate. 


Time (sec) 

4000 
3500 
3000 
2500 
2000 


1500 

1000 


500 


0 

2 plates Cantilever Cylinder Cylinder 2 Cylinders 

plates into plate into wall 

HH Cray case (2) Is Symult case (2) ■■ Symult case (1 ) 



Figure 8. Comparison of execution times for case(l) (improved parallel contact) 
and case (2) (no parallel contact) (Ref. 4). 
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and case (2) (no parallel contact) (Ref. 4). 
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DESIGN ISSUES 


Currently it takes from three to six months to develop a validated full vehicle finite element model 
that can be effectively used for design direction. With the current push to reduce the automotive vehicle 
design cycle into the 24-month time range, this model building cycle is clearly unacceptable. Thus, even if 
we cut the execution time for the codes to minutes rather than hours, we still may not have significantly 
affected the design process. Although fully automatic mesh generators for two-dimensional fields have 
been around for almost ten years (for both triangles and quads), they have not been effectively used for 
assembled structures. This suggests that more work needs to be focused on this critical problem. 


Design issues 


Solution time 
10-20 hrs 



Model building 
& validation 

3-6 months 


Figure 9. Significance of improved modeling methods. 
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DESIGN ISSUES 


Finally, although linear elastic structural optimization is one of the most developed fields of 
engineering design research, there have been a mere handful of papers in the area of crashworthiness 
optimization. Based on work in linear elastic structural optimization, there are two potentially fruitful areas 
of research. First, there is the question of efficient calculation of sensitivity information. . This is made 
doubly difficult because of the presence of nonlinearities introduced by contact. Second is the idea of 
constructing inexpensive, robust local design space approximations. In the example shown below, a 
lumped mass model was used as the local approximation and the finite element model was used only at a 
few number of points in the design space. 


Front lower rail 



Force - deformation 
curve for front lower 
rail 


Force 


Deformation 


Spring - mass model 



Figure 10. Optimum design with crashworthiness constraints (Ref. 5). 
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R AND G PARAMETERS: ATB SIMULATOR UNLOADING BEHAVIOR 


This figure demonstrates how unloading is handled in the ATB occupant simulator for contact 
planes and seatbelts. Two parameters, R and G, are specified. R is the ratio of rebound energy to initial 
energy and G is the ratio of permanent deformation to maximum deformation. The program unloads along 
a quadratic curve which intersects the x-axis at Dl, the point defined by G. The quadratic function is 
chosen to match the requested value of R as closely as possible while satisfying the G condition exactly. 
This puts an upper limit on R, which for linear loading and unloading is seen to be 1.0 G. Excessive 
rebound energy from the seatbelts prompted an investigation of the program s sensitivity to these 
parameters. 



G = Ratio of Permanent to 
Maximum Deformation 

= Dl / D2 

= 0.85 (in above example) 


R = Ratio of Rebound Energy 
to Initial Energy 

= [ 1.0 - G ] for Linear 
Loading and Unloading 


Figure 1 
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MOTION OF SPHERE RELATIVE TO VEHICLE AFTER VEHICLE DECELERATES 


The sketch at the top of this figure shows an ATB model developed by the consulting firm, David 
James, Ltd. It is simply a sphere, restrained by a single belt, and sliding on a frictionless, decelerating 
plane (or vehicle). From their study they concluded that the ATB R and G parameters were essentially 
nonfunctional for the seatbelt. They passed their results on to UVA where initial simulations seemed to 
confirm their conclusions. However, changing the integration step-size caused a significant change in the 
results, particularly in the rebound energy of the sphere. A systematic study of step-size sensitivity was 
undertaken. The sketch at the bottom shows an equivalent restraint situation with the belt replaced by a 
contact plane. This model was used for comparison to the belt model. 



Figure 2 
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SPHERE DATA FOR BELT RESTRAINT 


This figure displays the family of curves that were generated by lowering the integration step-size 
from 1.0 msec (the value used in the initial R and G investigation) to 0.01 msec. The G value was 0.85. 
As the 1 .0 msec curve shows, the final rebound velocity is indeed greater than the maximum forward 
velocity. The energy return for the smallest step-size, however, is approximately 0.16, which is near the 
maximum for linear loading and unloading. Furthermore, the curves do appear to be converging. 
Unfortunately, the step-size required for convergence is extremely small, much smaller than is normally 
necessary for motion of an unbelted occupant. The conclusion of this study is that the belt algorithm does 
work, including the R and G parameters, but very small step-sizes are needed to get reasonable results. 



— • DT=0.100 

— DT=0.050 

— DT=0.025 

— DT=0.010 


Figure 3 
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SPHERE DATA FOR WALL IMPACT 


To check that the integration problems were limited to the seatbelt, an equivalent model was set up 
as described m Fig. 2, with the belt replaced by a contact plane (or wall). This figure shows that the wall ’ 
impact results are nearly identical to the best belt results and that there is no discemable difference between 
the wall results for 1.0 msec and for 0.01 msec. In fact, the two cases exhibited four-digit agreement. 
Thus, the basic integration method is apparently sound and problems are limited to the belt algorithm. 



TIME (MSEC) 

x-x DT=1.000 
DT=0.01 0 
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Figure 4 



HEAD ACCELERATION VS. BELT SLACK 
X-DIRECTION (MEASURED) 


This figure shows head acceleration curves from two UVA impact sled tests which used the same 
parameters except for the amount of slack in the shoulder belt. Of particular interest is the secondary peak 
that occurs near 85 msec for the curve generated by the run with 3 inches of slack. This peak had been 
noticed before in sled tests and in simulations, but had been considered not to have any significance since 
its magnitude was always much lower than that of the primary peak. It appears, however, that the 
secondary peak may be associated with experimental parameters such as belt slack. Thus, while they are 
not of great importance from the standpoint of injury prediction, this and other secondary features have the 
potential to be of useful diagnostic value if their appearance in the pulse can consistently be associated with 
known parameters. 



Time (msec) 


Figure 5 
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HEAD ACCELERATION VS. BELT SLACK 
X-DIRECTION (SIMULATED) 


This figure shows ATB simulator results which are similar to the sled measurements. The shape 
and timing of the pulses need some improvement before they are compared directly to the sled results, but 
the appearance of a secondary peak resulting from the introduction of belt slack is clearly demonstrated. 




HEAD ACCELERATION VS. BELT SLACK 
X-DIRECTION (SIMULATED) 


This figure shows that the secondary peak is not an anomaly associated with a particular value of 
belt slack. It appears for higher and lower values as well, and demonstrates a predictable sensitivity to the 
amount of slack. Note that it shows greater sensitivity to the slack parameter than does the primary peak, 
particularly in terms of timing and the depth from the peak to the following trough. This increased 
sensitivity may be exploitable when diagnosing test results. 



Figure 7 


163 





/V 1 I 

N94-19474 


Current Capabilities for Simulating the Extreme 
Distortion of Thin Structures Subjected 
to Severe Impacts 


Samuel W. Key 
KEY Associates 
Albuquerque, NM 



PRECEDING PAGE BLANK NOT FILMED 


165 




SUMMARY 


The explicit transient dynamics technology in use today for simulating the impact and subsequent 
transient dynamic response of a structure has its origins in the "hydrocodes" dating back to the late 1 940’s. 
The growth in capability in explicit transient dynamics technology parallels the growth in speed and size of 
digital computers. Computer software for simulating the explicit transient dynamic response of a structure 
is characterized by algorithms that use a large number of small steps. 

In explicit transient dynamics software there is a significant emphasis on speed and simplicity. The 
finite element technology used to generate the spatial discretization of a structure is based on a compromise 
between completeness of the representation for the physical processes modelled and speed in execution. 
That is, since it is expected in every calculation that the deformation will be finite and die material will be 
strained beyond the elastic range, the geometry and the associated gradient operators must be 
reconstructed, as well as complex stress-strain models evaluated at every time step. As a result, finite 
elements derived for explicit transient dynamics software use the simplest and barest constructions 
possible for computational efficiency while retaining an essential representation of the physical behavior. 
The best example of this technology is the four-node bending quadrilateral derived by Belytschko, Lin and 
Tsay (1984). 

Today, the speed, memory capacity and availability of computer hardware allows a number of the 
previously used algorithms to be "improved." That is, it is possible with today's computing hardware to 
modify many of the standard algorithms to improve their representation of the physical process at the 
expense of added complexity and computational effort. 

The purpose of this presentation is to review a number of these algorithms and identify the 
improvements possible. In many instances, both the older, faster version of algorithm and the improved 
and somewhat slower version of the algorithm are found implemented together in software. Specifically, 
the following seven algorithmic items are examined: 

1) The invariant time derivatives of stress used in material models expressed in rate form. 

2) Incremental objectivity and strain used in the numerical integration of the material models. 

3) The use of 1 -point element integration versus mean quadrature. 

4) Shell elements used to represent the behavior of thin structural components. 

5) Beam elements based on stress-resultant plasticity versus cross-section integration. 

6) The fidelity of elastic-plastic material models in their representation of ductile metals. 

7) The use of Courant subcycling to reduce computational effort. 


PRBQOWG PAGE BLANK NOT FH-ME 




167 



INVARIANT TIME DERIVATIVE OF STRESS 


To account for the fact that bodies subjected to large displacements undergo significant rigid body 
rotations, material models in rate form rely on an invariant time derivative of the stress. For example, a 
linear elastic material expressed in terms of rates has the form: 


= C rsmn d , 


( 1 ) 


where cojbn is a rotation rate, dmn is the stretching and C rsmn y the material's elastic modulus. To date the 
majority of simulations have used the Jaumann invariant derivative which uses the spin for the rotation 


rate, co*™ = ('Ok.m - ^W)/2. 

This formulation is very sensitive to shear strains that are greater than 20% in the presence of rotations 

(Dienes, 1979). . . . , . - _ _ 

A more accurate invariant time derivative is the Green-Mclnnis derivative which is based on the 


rotation from the polar decomposition of the deformation gradient, = R \ n R : , where F™ = Rmn u nc 
(Dienes, 1979). 

Computing the uniform shearing of a 1 x 1 coupon demonstrates the difference in predicted behavior 
using kinematic hardening plasticity for these two invariant time derivatives (Fig. 1). 
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(Note: In (£/t 0 ) = 5 => C = 10) 

Figure 1 - Uniform shearing of a 1 x 1 coupon 



INVARIANT TIME DERIVATIVE OF STRESS (Cont'd.) 


An examination of the shear stress as a function of shear strain produces very interesting results (Fig. 
2). The results based on using the Green-Mclnnis derivative in Fig. 2 are denoted by "Dienes." As can be 
seen, extreme distortions using the Jaumann derivative lead to physically unrealistic stress variations that 
change sign. The monotone increasing curve denoted by "Dienes" exhibits a physically realistic stress- 
strain representation. This anomalous behavior exhibited by the Jaumann derivative must be avoided in a 
simulation if practical results are to be obtained. 

For the sake of reliability, today's software will offer the Green-Mclnnis form of the invariant time 
derivative of stress along side of the Jaumann derivative for material models based on rate formulations. 



Figure 2 - Kinematic hardening plasticity predictions of shear 
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INCREMENTAL OBJECTIVITY FOR STRAIN 


In numerical simulations differential expressions are replaced with incremental expressions. 
Incremental expressions for such things as rotation and strain increments can produce errors even for small 
steps if the properties of the differential expression are lost. For example, when the time derivative of an 
orthogonal rotation is approximated, it is very easy to lose the orthogonal properties associated with the 
differential form. 

The same thing can happen when integrals of strain rates are replaced by sums of strain increments, 
particularly when finite strains are involved. Traditionally, a strain increment is obtained from the 

Stretching with A Emn ~ At d mn > d mn — (y)m,n + V n,m )/2* 

This approach to generating strain increments can be sensitive to cyclic strain histories in the presence 
of rotations. 

A more accurate, but computationally more costly approach, is to base the strain increment on the 
symmetric stretch tensor U obtained from the polar decomposition of the deformation gradient, 

pm - ftmn jj This latter approach has been termed strong objectivity by M. Rashid (1992), who 

(X fiCC 1 1 

refers to the traditional approach as weak objectivity. 

Computing the hoop vibration of a rotating ring provides a good example of the contrasting results that 
can come from the use of the faster, classical formulations that are weakly objective and the more accurate 
and costly formulations that are strongly objective. Figure 3 shows an elastic ring given an initial angular 
velocity of 4,000 radians per second. 



Figure 3 - An elastic ring with an initial 4,000 rad/s angular velocity 
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INCREMENTAL OBJECTIVITY FOR STRAIN (Cont'd.) 


The ring rotates approximately 360 degrees in 1 .5 milliseconds. During that time the ring oscillates 
radially five times (Fig. 4). Of particular interest is the increase in kinetic energy with time (Fig. 4). This 
is a closed system to which no additional energy is being added. The increase in kinetic energy reflects the 
accumulation of errors from the strain increments in the classical formulation that is weakly objective in 
spite of the very small steps in the algorithm. 


Rotating Ring (Omega - 4,000 rads/sec) HXEL Oaumann, Elastic) FMA-3D V07. 32 31-3UL-1992 09:50:09 
Strain Based on E*n - (t(n*l)-t(n)) x <Vm,n * Vn,m>, WeaV Objectivity 



Figure 4 - Kinetic energy response with only weak objectivity 
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INCREMENTAL OBJECTIVITY FOR STRAIN (Cont'd.) 


The same computation using a formulation that is strongly objective in the sense of Rashid shows the 
expected response (Fig. 5); namely, a kinetic energy that oscillates between two limits that are constant in 
time. This latter calculation shows an exemplary energy exchange between kinetic energy and internal 
strain energy as the ring rotates and oscillates radially. 

In spite of the increased expense, today's software should offer as the default a strain increment 
formulation meeting the requirements for strong objectivity. 


Rotating Ring (Onega - 4,000 rads/ sec) HXEL (Polar Deconp, Elastic) l-JUL-1992 08:51:45 
Strain Based on N, Rashid's Definition of Strong Objectivity 



Figure 5 - Kinetic energy response with strong objectivity 
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ELEMENT INTEGRATION 


In the interest of speed, explicit transient dynamic computer programs routinely use constant stress 
elements. In effect, only the uniform strain states of the element are used to represent the structures 
behavior. Strain states and, consequently, stress states that vary over the element are ignored. Such 
approximations result in complex, short wavelength, stress-free deformation patterns that are called 
hourglass modes. 

The quickest way to obtain a constant stress element is to evaluate the integrals over the area or volume 
of the element defining the gradient and divergence operators with a single integration point; hence, the 
terminology 1-point integrated elements. Unfortunately, such elements can fail the Iron s patch test. That 
is, a collection of irregularly-shaped elements subjected to a linear motion on the boundary will not 
produce a constant state of strain in the interior and, therefore, the collection of irregularly-shaped elements 
will not produce a constant state of stress. 

A more accurate approach is to evaluate the integrals over the area or volume of the element defining 
the gradient and divergence operators exactly for a constant state of stress, effectively a projection 
calculation. The result is a much more reliable and well-behaved simulation, albeit requiring the execution 
of a greater number of algebra expressions. The mean quadrature elements while still constant strain 
elements, pass the Iron's patch test. 

The greatly reduced excitation of the hourglass modes and the greatly improved hourglass control 
offered by the mean quadrature integration over the 1 -point integration virtually renders obsolete the older 
1-point integration technology. (Tn the special case of two-dimensional continuum elements, both 
approaches yield the same results.) 

Current users of explicit transient dynamic software should only expect to use 1 -point integrated 
elements in place of elements based on mean quadrature to obtain "quick and dirty results. 
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QUADRILATERAL SHELL ELEMENTS 


In explicit transient dynamics software there is a significant emphasis on speed and simplicity. That is, 
since it is expected in every calculation that the deformation will be finite and the material will be strained 
beyond the elastic range, the geometry and the associated gradient operators must be reconstructed, as well 
as, complex stress-strain models evaluated at every time step. As a result, finite elements derived for 
explicit transient dynamics software use the simplest and barest constructions possible for computational 
efficiency while retaining an essential representation of the physical behavior. The best example of this 
technology is the four-node bending quadrilateral derived by Belytschko, Lin and Tsay (1984). 

The BLT element is based on a constant stress assumption coupled to a flat plate geometric 
approximation of what is usually a warped geometry (the element's geometry is warped when the four 
nodal points do not define a flat surface). In certain situations the BLT element exhibits shortcomings in 
representing the response of a structure. Two examples are the bending of a twisted beam (Example 1), 
and the response of a hemisphere loaded by opposing forces across the hemisphere's diameter (Example 
2 )- 

A four-node bending quadrilateral has been developed that exhibits the expected behavior in these two 
examples. The element has two properties that provide the expected response: 

1) a warping deformation that possesses proper structural stiffness as opposed to being an hourglass 
mode, and 

2) a derivation that is based on the actual geometry of the element as opposed to treating the geometry 
as flat. 

Example 1 . The tip-loaded, twisted cantilever beam is an example of a structure with a nonplanar 
geometry that can be nontrivial to reproduce with the finite element technology available today in explicit 
transient dynamic software (Fig. 6). 



Figure 6 - Twisted beam with a constant tip load 
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QUADRILATERAL SHELL ELEMENTS (Cont'd.) 


The tip load is applied suddenly at time equal to zero and held constant in magnitude and orientation. 
Based on a length of 12.0 in., a width of 1.1 in., a thickness of 0.32 in. and the selected properties (a 

Young's modulus of 2.9 x 10 7 psi, a Poisson's ratio of 0.22, and a density of 2.5 x 10 4 lbf-sec 2 /in 4 ), 
the static deflection equals 0.005424 in. (MacNeal and Harder, 1985). The fundamental period equals 8.0 
milliseconds. 

As can be seen from Fig. 7, the BLT element with zero hourglass stiffness predicts an unacceptably 
large deflection amplitude and response period; see also the results reported by Belytschko, Wong and 
Chiang (1989). 


Bea« with a 90 Degree Spiral Twist (BelytschKo-Tsay Shell) 19:52:22 17-JUL-91 



Figure 7 - Beam tip motion prediction using Belytschko-Tsay shell 
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QUADRILATERAL SHELL ELEMENTS (Cont'd.) 


Example 2. A hemisphere loaded by two sets of diametrically opposed forces in the plane of the 
equator, separated by 90 degrees and alternating in sign, is a problem in which both bending and 
membrane deformation occur. The loads enter the structure by generating moments including warping or 
twisting moments in the elements adjacent to the loads. 

The loads are applied suddenly at time equal to zero and held constant in magnitude and orientation 
(Fig- 9). 



Figure 9 - Quadrant of a hemisphere loaded across two diameters 
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QUADRILATERAL SHELL ELEMENTS (Cont’d.) 


Based on the hemisphere's radius of 10.0 in., a thickness of 0.04 in. and the selected properties (a 

Young's modulus of 6.825 x 10 7 psi, a Poisson's ratio of 0.3, and a density of 2.5 x 10 4 lbf-sec 2 /in 4 ), 
the static deflection equals 0.094 inches (MacNeal and Harder, 1985). 

As can be seen from Fig. 10, the BLT element with zero hourglass stiffness predicts an unacceptable 
cyclic accumulation of displacement; see also the results reported by Belytschko, Wong and Chiang 
(1989). 


Diameter Forces On A Hemisphere (Belytschko-Tsay Shell) 19:56:34 17-JUL-91 



Figure 10 - Radial motion prediction using Belytschko-Tsay shell 
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QUADRILATERAL SHELL ELEMENTS (Cont’d.) 


The technology developed by KEY Associates predicts a steady periodic (constant amplitude) result 
that is within 0.2% (Maximum_Internal_Energy / 4.0) of the expected deflection (0.094 inches) (Fig. 1 1 ). 

Remarks. This element is suitable for both large in-plane and bending strains. These particular 
calculations were carried out without any hourglass control suggesting this element will not be overly 
sensitive to the development of hourglass distortions. In addition, this technology is based on six (6) 
degrees of freedom per nodal point meaning that no further numerical constraints or adaptations are needed 
to eliminate the drilling" rotation to obtain satisfactory results. Modeling such additional physical features 
as edge beams, or modeling folded plate structures does not present any particular difficulty. 

The low level of loading in both of these examples results effectively in infinitesimal strains. The 
material response is linear elastic. However, both of these examples require a proper computation of the 
gradient and divergence operators to obtain the correct results. The examples are sensitive indicators of the 
correctness of the representation for the element geometry and the element twisting stiffness. 

While the BLT four-node bending quadrilateral remains a computationally efficient element, there are 
numerous applications for which an accurate representation of the warped geometry and the twisting 
stiffness is essential to obtaining a satisfactory result. Without a doubt, up-to-date software should contain 
an efficient four-node bending quadrilateral with the capabilities discussed here. 


Diameter Forces On A Hemisphere (Pl«tel6, 10/2/1) JAWS -3D V06. 01 16 -NOV -1991 08: 32: 12 



Figure 1 1 - Radial motion prediction (= I E/4) using improved four-node shell 
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QUADRILATERAL SHELL ELEMENTS (Cont'd.) 


The technology developed by KEY Associates predicts a steady periodic (constant amplitude) result 
that is within 0.2% (Maximum_Internal_Energy / 4.0) of the expected deflection (0.094 inches) (Fig. 11). 

Remarks. This element is suitable for both large in-plane and bending strains. These particular 
calculations were carried out without any hourglass control suggesting this element will not be overly 
sensitive to the development of hourglass distortions. In addition, this technology is based on six (6) 
degrees of freedom per nodal point meaning that no further numerical constraints or adaptations are needed 
to eliminate the "drilling" rotation to obtain satisfactory results. Modeling such additional physical features 
as edge beams, or modeling folded plate structures does not present any particular difficulty. 

The low level of loading in both of these examples results effectively in infinitesimal strains. The 
material response is linear elastic. However, both of these examples require a proper computation of the 
gradient and divergence operators to obtain the correct results. The examples are sensitive indicators of the 
correctness of the representation for the element geometry and the element twisting stiffness. 

While the BLT four-node bending quadrilateral remains a computationally efficient element, there are 
numerous applications for which an accurate representation of the warped geometry and the twisting 
stiffness is essential to obtaining a satisfactory result. Without a doubt, up-to-date software should contain 
an efficient four-node bending quadrilateral with the capabilities discussed here. 


Diameter Forces On A Hemisphere (Pl«tel6, 10/2/1) JAWS-3D V06 01 16-N0V-1991 08:32:1 2 



Figure 1 1 - Radial motion prediction (= I E/4) using improved four-node shell 
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BEAM ELEMENTS 


Beam elements, while on the surface appearing to be very elementary structural components, can be 
very intense from a computational standpoint when finite strains and inelastic material behavior are 
modeled. A portion of the computational complexity comes from the need to evaluate the stress-strain 
behavior at many points over the cross section (16 integration stations is a typical number). Membrane, 
bending and torsional stress resultants are produced from weighted sums over these points. 

A very fast and direct method of obtaining stress resultants for elastic-plastic material behavior is to 
introduce stress-resultant plasticity. In effect, the beam cross-section is either completely elastic or 
completely plastic. For simulations in which the beam is deformed well into the plastic range, 
considerable computational effort can be eliminated with stress-resultant plasticity (the evaluation of 
elastic-plastic stress-strain models at numerous individual points over the cross section is not required). 

If only a portion of the beam cross section yields plastically, stress-resultants plasticity will not capture 
any of the detail. 

Quality software seeking to provide options for both rapid results and accurate results will offer both 
beam formulations. 



COURANT SUBCYCLING 


The central difference time integrator is only conditionally stable. That is, the integration procedure 
must be used with a time step that is less than the critical time step. The critical time step is very nearly 
equal to the minimum transient time for a disturbance to cross between any two nodal points in the mesh. 
The algorithm is very effective for an excitation that uses the maximum resolution the mesh can provide, 

for example, a shock wave. .... . . . , , „ , „ 

There are two very common circumstances when an explicit time integration algorithm becomes less 

efficient: first, when there are significant differences in the spacing between nodal points over the mesh 
and, second, when there are significant differences in material stiffness over the mesh. It is not 
uncommon to have differences in critical time steps as much as a hundred to one over a mesh, l he 
consequence is that the central difference time integration must function with the smallest critical time step. 

Fortunately, Courant subcycling may be introduced in order to recover much of uie efficiency ottered 
bv explicit algorithms. In Courant subcycling each finite element is integrated with the largest time step 
permitted by the local critical time step. Thus, small, stiff elements are integrated with many small time 
steps, and elsewhere in the mesh large, soft elements are integrated periodically with their larger time 

^The benefit is significant since the time to perform a calculation drops and the answers are more 
accurate since the central difference algorithm performs as close to the critical step as is practtcab e 
everywhere. Figure 12 shows an example of a domain in which elements are integrated with two different 
time steps. In this case the ratio is only 2 to 1, but the principle is demonstrated. 

Clearly, software intended for a wide variety of applications should provide Courant subcycling. 



Figure 12 - Element integration bins for Courant subcycling 
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PLASTICITY IMPLEMENTATIONS 


Traditional elastic-plastic material models are based on a linear representation of plastic hardening and 
radial-return numerical integration. While an exceptionally effective approach for large plastic strains, it is 
somewhat approximate for moderate plastic strains (plastic strains less than a ten times elastic strains. 
Ductile metals can exhibit a very rounded" stress-strain curve when first yielding plastically. An elastic- 
linear strain hardening model with its transition from elastic to plastic response at a sharp comer does not 
mimic a gradual transition to large plastic strains. 

One should expect a variety of material models in a crash simulation program including the basic linear 
strain hardening model and a model for plastic yielding that has a smooth transition from elastic to plastic 
response. For small plastic strains the models available should include sub-stepping where the time step 
within the material model is broken into smaller steps to control the numerical integration error. 
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SOFTWARE SURVEY 


The following two tables (Tables 1 and 2) are a modest survey of the crash simulation software 
represented by the speakers at this NASA Workshop on Crashworthiness Technology. In brief, the 
software surveyed is as follows: 

1. FMA-3D is a full-feature program available commercially from KEY Associates, 1851 Tramway 
Terrace Loop NE, Albuquerque, NM 87122 that include subcycling. 

2. DYNA-3D is a full-featured program available through a technology sharing agreement from the 
Lawrence Livermore National Laboratory, Attn: Robert Whirley, P.O. Box 808, Livermore, CA 94550. 

3. PRONTO-3D is a limited-feature program available with a license agreement from Sandia National 
Laboratories, Attn: Steven Attaway, Organization 1425, Albuquerque, NM 87185. 

4. SUPERWHAMS is a full-featured program soon to be available commercially from KBS2, Suite 
1 12, 455 South Frontage Road, Burr Ridge, IL 60521. 

5. DYCAST is a limited-feature program based on an implicit integration algorithm available 
commercially from Grumman Aerospace Corporation, Attn: Allan Pifko, MS A08-35, Bethpage, NY 
11714. 


Table 1. Technology Availability 



FMA-3D 

DYNA 

PRONTO 

Time Derivative 

Jaumann 

Yes 

Yes 

No 

Green-Mclnnis 

Yes 

No 

Yes 

Strain Increment 

Weak Objectivity 

Yes 

Yes 

No 

Strong Objectivity 

Yes 

No 

Yes 

Element Integration 

1-Point 

No 

Yes 

No 

Mean Quadrature 

Yes 

Yes 

Yes 

4-Node Shells 

BLT Quad 

Yes 

Yes 

Yes 

Improved Quad 

KEY 

EWQ 

No 

Beam Elements 

Illyshin Plasticity 

(i/0 

Yes 

No 

Integrated 

Yes 

Yes 

NO 

Plasticity 
Linear Hardening 

Yes 

Yes 

Yes 

w/ Radial Return 

Smooth Hardening 

Yes 

"Yes” 

No 

w/ Sub-Stepping 

Time Subcycling 

Yes 

No 

No 


Table 2. Technology Availability 



SUPER WHAMS 

DYCAST 

Time Derivative 



Jaumann 

Yes 

(no solid 

Green-Mclnnis 

No 

elements) 

Strain Increment 



Weak Objectivity 

Yes 

(no solid 

Strong Objectivity 

No 

elements) 

Element Integration 



1-Point 

No 

(no soiid 

Mean Quadrature 

Yes 

elements) 

4-Node Shells 



BLT Quad 

Yes 

(A l/d&r) 

Improved Quad 

BWCQ 

(inf strain) 

Beam Elements 
Illyshin Plasticity 

No 

No 

Integrated 

Yes 

Yes 

Plasticity 



Linear Hardening 
w/ Radial Return 

Yes 

No 

Smooth Hardening 
w/ Sub-Stepping 

No 

Yes 

Time Subcycling 

Yes 

Implicit 
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CRASHDYNAMICS WITH DYNA3D: CAPABILITIES 

AND RESEARCH DIRECTIONS 


This paper discusses the application of the explicit nonlinear finite element analysis code DYNA3D to 
crashworthiness problems. Emphasized in the fu st part of this work are the most important capabilities of 
an explicit code for crashworthiness analyses. The remainder of the paper discusses the areas with 
significant research promise for the computational simulation of crash events. 
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WHAT IS CRASHWORTHINESS? 


In the present context, crashworthiness is defined as the study of vehicle survivability under the impact 
with another object. It is the focus on survivability which differentiates a crashworthiness analysis from 
other kinds of vehicle dynamic analysis. Crashworthiness analyses include the simulation of multi-car 
collisions, impacts of automobiles into highway barriers, as well as aircraft collisions with terrain or 
airborne objects such as birds. 

The Laboratory's involvement with crashworthiness began in a passive role in the mid 1980's, when 
many carmakers began using our DYNA3D software for crashworthiness analysis and providing us 
feedback on its performance. Recently, the Laboratory has taken a more active role in crashworthiness 
analysis through a DOT highway project, the DOE Auto Initiative, and involvement in the Electric Vehicle 
Consortium. 


Crashworthiness is the study of vehicle survivability under impact 
with another object. 

• Auto/Auto 

• Auto/Highway barrier 

• Aircraft/Terrain 

• Aircraft/Bird 

LLNL's involvement in crashworthiness began with DYNA3D, and 
more recently: 

• DOT Highway Project 

• DOE Auto Initiative 

• Electric Vehicle Consortium 
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CHARACTERIZATION OF CRASH EVENTS 


Crash events are characterized by highly dynamic transient structural response of 5-100 millisecond 
duration. Mechanical contact conditions are almost always present. In addition, vehicle crash events 
frequently involve large deformations of thin structures and highly inelastic material response. 1 hese 
characteristics suggest explicit nonlinear finite element analysis as an effective tool for crashwort mess 
analysis. 


Crash events usually involve: 

• Transient dynamic effects (5- 1 00 ms duration) 

• Mechanical contact conditions 

• Large deformations of thin structures 

• Inelastic material response 

These characteristics often motivate the use of an explicit nonlinear 
finite element code for crashworthiness analysis. 
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CRASH STUDIES USE FINITE ELEMENT MODELS AT SEVERAL LEVELS 


Although many people identify crashworthiness analysis with very large full-body models, effective 
crash analysis is often conducted with a hierarchy of models. These models range from single 
components to subassemblies containing a small number of components to full vehicle models. The 
appropriate level of model complexity is selected based on considerations of simulation accuracy versus 
the time required to construct the model, run the analysis, and interpret the results. 


Component 



rail collapse 


Subassembly 


Full Vehicle 




B 1 -B side impact frontal impact 


Tradeoff is simulation accuracy vs. time to: 

• Construct model 

• Run analysis 

• Evaluate results 
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KEY DYNA3D CAPABILITIES FOR CRASHWORTHINESS MODELING 


The following viewgraphs discuss some of the capabilities of DYNA3D which have been found 
important for crashworthiness modeling. The selection of these features as key capabilities was based on 
both our experience at LLNL and on extensive feedback received from the DYNA3D outside user 
community. Viewed collectively, this list of capabilities depicts a general computational simulation 
capability which is representative of the current state-of-the-art in crashworthiness software used for 
production analysis, including both DYNA3D and its derivatives as well as other explicit FE crash codes. 


Based on our experience + feedback from user community: 

• Rigid materials/bodies 

• Good structural elements (shells and beams) 

• Fast unilateral contact 

• Robust general contact algorithms 

• Variety of constitutive models 

• Failure modeling 
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KEY CAPABILITIES - GOOD STRUCTURAL ELEMENTS 


Most large crashworthiness models use some beam elements to represent vehicle structural members. 
DYNA3D contains the Hughes-Liu beam element and the Belytschko-Schwer beam element. The 
Hughes-Liu beam uses numerical integration in the cross section and contains a moveable reference 
surface, and is therefore more general but also more expensive. More commonly used in crash problems 
is the resultant-based Belytschko-Schwer beam element. 


Beams are also frequently useful in vehicle modeling. 

• Beam elements: 

• Hughes-Liu 

• Belytschko-Schwer 

• Simple truss 

• One point integration along the length 

• Hughes-Liu uses variable-order numerical integration 
in the cross-section for accurate nonlinear material 
behavior 

• Belytschko-Schwer is simple resultant formulation 
(much less expensive) 
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KEY CAPABILITIES - RIGID MATERIALS/BODIES 


The capability of modeling portions of a vehicle as rigid and other portions as deformable in the same 
analysis is quite useful in crashworthiness analysis. So-called "rigid materials" can be used to correctly 
represent the mass and inertial properties of part of a vehicle which experiences little deformation at a small 
fraction of the cost required if all deformable materials were used. Another common application of rigid 
materials is the inexpensive definition of a complex curved surface for contact calculations, such as a rigid 
pole or partial barrier to be struck by an oncoming vehicle. Finally, rigid materials serve as a powerful 
model debugging feature. Portions of the model can be easily changed from deformable to rigid, thereby 
isolating regions of the model which may be the source of difficulties. 


• Define all elements of a specified material as composing a rigid 
body (inertial properties automatically computed from geometry 
or specified by analyst) 

• Greatly reduces cost compared to deformable elements 

• Useful for: 

• representing mass properties for parts of vehicle 
experiencing little deformation 

• defining a complex rigid surface for contact 

• model debugging 



KEY CAPABILITIES - GOOD STRUCTURAL ELEMENTS 


A variety of good structural elements are essential for vehicle crashworthiness analysis. DYNA3D 
contains a library of quadrilateral and triangular shell elements, each with unique capabilities. All elements 
share the common features of one-point in-plane integration with stabilization of zero-energy modes, and 
variable order numerical integration in the through-thickness direction. 


A family of shell elements with common features gives'versatility. 

• One point in-plane integration with stabilization 

• Variable-order numerical integration thru-thickness 

• Shell elements (4-node quads and 3-node triangles): 

• Hughes-Liu 

• Belytschko-Tsay 

• C° triangle 

• BCIZ triangle 

• Membrane (derived from B-T) 

• YASE 
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KEY CAPABILITIES - GOOD STRUCTURAL ELEMENTS 


Most large crashworthiness models use some beam elements to represent vehicle structural members. 
DYNA3D contains the Hughes-Liu beam element and the Belytschko-Schwer beam element. The 
Hughes-Liu beam uses numerical integration in the cross section and contains a moveable reference 
surface, and is therefore more general but also more expensive. More commonly used in crash problems 
is the resultant-based Belytschko-Schwer beam element. 


Beams are also frequently useful in vehicle modeling. 

• Beam elements: 

• Hughes-Liu 

• Belytschko-Schwer 

• Simple truss 

• One point integration along the length 

• Hughes-Liu uses variable-order numerical integration 
in the cross-section for accurate nonlinear material 
behavior 

• Belytschko-Schwer is simple resultant formulation 
(much less expensive) 
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KEY CAPABILITIES - UNILATERAL CONTACT (RIGID WALLS) 


Unilateral contact, known as "rigid walls" in DYNA3D, is another widely used feature for 
crashworthiness analysis. This option allows a simple definition of a rigid plane for contact, and is often 
used to simulate vehicle barrier crash tests. Unilateral contact offers execution speed and modeling 
simplicity as advantages over discretizing the barrier and using general two-surface contact to treat the 
impact conditions. 


• Applicable for deformable vehicle impact onto rigid 
barrier 

• Much less expensive than discretizing target and using 
general contact formulation 

• Often used for auto frontal crash 
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KEY CAPABILITIES - ROBUST GENERAL CONTACT 


Robust general contact algorithms may be the most important capability in a crashworthiness code. 
DYNA3D contact is based on a slave node on master segment concept, and a two-way symmetric 
treatment is used to eliminate any bias in the calculations. Important components of this capability are the 
treatment of contact between solid and shell elements, between beams and other element types, and single 
surface (self) contact. 


• Slave node on master segment, symmetric treatment, 
penalty-method based 

• Two-surface algorithm: 

• arbitrary interaction of rigid and deformable bodies 

• general treatment of solids and shells, beams by 
node only 

• incremental search 

• Single-surface algorithm (self-contact) - buckling and 
folding 
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KEY CAPABILITIES - FAILURE MODELING 


Failure modeling is evolving as a useful tool for crashworthiness modeling. Current approaches in 
DYNA3D include nodal constraint release options, material-based failure, and element deletion based on a 
failure criterion. Although these approaches have repeatedly proved their value to skilled analysts, they are 
somewhat ad-hoc and lack theoretical basis in general settings. 


• Nodal constraint release 

• "tie-breaking shell slidelines" 

• "tied node sets with failure" 

• Material-based failure - element no longer carries stress 
(deviatoric or total) 

• These approaches are somewhat ad-hoc and may be mesh- 
dependent, but have proven useful to skilled analysts 

• Significant improvements are needed here. 
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CRASHWORTHINESS SIMULATION 
ACTIVITIES AT LLNL 


Crashworthiness efforts at LLNL include both computational methods research and applications. 


• Computational methods 

• YASE shell - resistant to hourglassing 

• SAND - adaptive slidesurfaces for failure modeling 

• Constitutive models and failure criteria - metals and 
composites 

• ParaDyn - massively parallel DYNA3D 

• Simulation techniques 



YASE SHELL ELEMENT 


The YASE shell element is a recent outgrowth of research at LLNL. This element is a four-node 
quadrilateral which extends ideas and the coordinate system of the Belytschko-Tsay element. In the YASE 
shell, the stabilization evolves directly from the formulation and there are no free parameters which must 
be chosen by the user. This element yields improvements in accuracy over the Belytschko-Tsay element, 
and is comparable in performance to more recent elements developed by Belytschko and coworkers. The 
YASE shell as implemented in DYNA3D is within 10% of the speed of the Belytschko-Tsay element and 
therefore does not represent a substantial increase in analysis cost. 


• Builds on ideas and coordinate system of Belytschko-Tsay 

• Four-node quad, resultant-based, element normal 

• Stabilization evolves directly from the formulation - no 
"tunable parameters" 

• Good coarse mesh accuracy, especially for in-plane 
bending 

• Speed competitive with Belytschko-Tsay (the fastest 
DYNA3D quad shell) 
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SAND-SLIDESURFACES WITH ADAPTIVE NEW DEFINITIONS 


A recently developed DYNA3D capability for modeling material failure on contact surfaces is SAND 
(Slidesurfaces with Adaptive New Definitions). Engineering analysts have found SAND useful in 
modeling bird impacts on windshields and airframes as well as vehicle impacts on soft "frangible" 
roadside structures. SAND acts as an algorithmic vehicle to implement failure in an explicit finite element 
code; the difficult task of defining appropriate failure criteria for different problem classes remains a 
significant area of research. 


• Models material failure on contact surfaces 

• Failed elements are deleted and contact surface adapts to 
new outer material boundary 

• Permits structural modeling using solid and shell elements 

• Allows penetration and failure modeling in Lagrangian 
framework 

• Useful for bird strikes on airframe, vehicle impacts on 
multiple soft barriers 



PARADYN - MASSIVELY PARALLEL DYNA3D 


The development of massively parallel explicit nonlinear finite element software is the goal of the 
LLNL ParaDyn Project. The move to massively parallel computing hardware will allow the solution of 
much larger problems with reduced turnaround times. It is important to realize, however, that MPP may 
not offer large speedups for today's small and moderate-sized problems due to difficulties with load 
balancing and multiprocessor overhead. The true potential of MPP-based analysis lies in the ability to 
incorporate levels of modeling detail which are totally infeasible with current computing hardware, and in 
the promise of improved price/performance across a range of machines. 


Major Objectives: 

• Develop a new generation of DYNA codes to take full 
advantage of teraflop computing platforms 

• Solve significantly larger problems in times 
commensurate with design cycles 

• Determine the architecture and parallel programming 
models best suited to explicit finite element methods 

• Architectures: CM-5, KSR, Intel, Cray MPPO 

• Programming Models: F90, HPF, Message Passing 

• Establish a community of university and industrial 
collaborators 
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DYNA3D IS USED TO MODEL THE UCD/CALTRANS TEST BOGEY. PROPER 
MATERIAL MODEL IS NEEDED FOR AGREEMENT WITH TEST DATA 


One recent LLNL crashworthiness study involved modeling a frontal crush unit tested on the UC 
Davis/Caltrans test bogey. The graph shows that agreement between the force-displacement curves 
measured in the experiment and predicted by DYNA3D was significantly improved when the forming-limit 
failure criterion was added into the DYNA3D model. Other LLNL crash activities include modeling the 
impact of the vehicles into various roadside barriers and objects such as light poles. 



"EI-Pr.Elastlc-PUstic: too (tiff 

"FLD"=FLDP Failure: Agrees with teat. 



WE ARE WORKING WITH FHWA AND 
U. ALASKA TO DEVELOP COMPLETE 
VEHICLE- BARRIER SIMULATIONS 
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DATA WAS SUPPLIED BY NHTSA FOR A 30 MPH FRONTAL IMPACT. THAT 
EVENING, WE RAN OUR SATURN MODEL "AS-IS" IN POST-PREDICTIVE MODE 


Another recent crashworthiness effort at LLNL consisted of taking our simplified Saturn model 
(developed in collaboration with J. Wekezer and shown on the lower right of this figure) and simulating a 
30 mph frontal crash into a rigid barrier. Test data was supplied by the National Highway Traffic Safety 
Administration (NHTSA) at three locations: front brake caliper, engine, and rear seat. TJis was not 
intended as a detailed modeling effort, but rather an appraisal of the usefulness of our highly simplified car 
model for representing gross response behavior in an auto crash scenario. 


Accelerations were supplied 
by NHTSA 

At three locations: 

Wheel 
- Engine 

* Rear Seat 

These were compared to the 
first and only DYNA3D ran 
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DYNA3D ACCELERATIONS OF ENGINE AND REAR SEAT AGREE WITH 
TEST DATA, ALLOWING FOR BUMPER CRUSH 


The simplified Saturn model was run once in DYNA3D for the 30 mph frontal crash, and the results 
IS2 npai ^ Wth te u St d £ to for j h . e en 8 ine block and rear seat locations. Reasonably good agreement was 
C -° nS,d j nng th , at the mode p. re ' dated the test data and was not tuned, and considering that the model 
was of a size and complexity to permit overnight analysis on an engineering workstation. Although clearly 

oSe^fh 1 f0r f d n t ^ ed r ful ^ b0dy , CI ? Shdyn 1 ami cs models ’ ±is does illustrate that useful results can Y 
be obtained with carefully defined simple finite element models. 


At Engine Block: 

- DYNA curve affect 30ms 
for hamper crush 

* DYNA shows less 
rebound due to rigid 
motor mounts 



At Rear Seat: 

- DYNA shows early peak 
due to rigid model aft of 
firewall 

- Overall Impulse and 
duration agrees, allowing 
for bumper crush 
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FUTURE DIRECTIONS 


Although substantial progress has been made in the computational modeling of vehicle 
crashworthiness in recent years, clearly much remains to be done. The computational modeling of 
composite materials, particularly under severe dynamic axial loads, remains an area in need of additional 
research. The general area of failure modeling, including both metals and composites, needs a m ore 
fundamental basis and more robust computational methods. The treatment of small scale structural details 
such as fasteners and spotwelds in a large crashworthiness model merits examination. Asm the case ot 
failure models, ad-hoc methods exist which usually work in the hands of skilled analysts, but lac 
theoretical foundation and are difficult to apply in new situations without substantial reliance on small-scale 

experiments. 

In conclusion, there is an opportunity for shared benchmark problems to be posed as a vehicle for 
communication among the research, code development, and engineering analysis communities in both the 
automotive and aerospace industries. These benchmarks could serve as illustrations of current needs in 
industry and as concrete objectives for researchers and code developers. These should be viewed as 
learning exercises and growth opportunities by all parties, however, and should not be viewed as software 

evaluation criteria. 


• Improved constitutive models for composite structural 
failure, especially axial compression - need increasing 
with growing use of composites in aircraft and vehicles 

• Improved progressive failure modeling - current 
procedures are ad-hoc but usually work. Need better 
basis and more robust methods 

• Special elements for modeling spotwelds and other small- 
scale features. Current tied-node capability (DYNA3D) 
is not sufficiently general and may excite hourglassing. 

• There is a strong need in both aerospace and automotive 
areas for shared benchmark problems of current relevance. 
These could be a strong vehicle for communication 
between crashworthiness analysis and research 
communities. 
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INTRODUCTION 


This report will present a brief overview of the transient dynamics capabilities at Sandia National 
Laboratories, with an emphasis on recent new developments and current research. In addition, the Sandia 
National Laboratories (SNL) Engineering Analysis Code Access System (SEACAS), which is a collection 
of structural and thermal codes and utilities used by analysts at SNL, will be described. The SEACAS 
system includes pre- and post-processing codes, analysis codes, database translation codes, support 
libraries, Unix shell scripts for execution, and an installation system. 

SEACAS is used at SNL on a daily basis as a production, research and development system for the 
engineering analysts and code developers. Over the past year, approximately 190 days of CPU time have 
been used by SEACAS codes on jobs running from a few seconds up to two and one-half days of CPU 
time. SEACAS is running on several different systems at SNL including Cray Unicos, Hewlett Packard 
PH-UX, Digital Equipment Ultrix, and Sun SunOS. 

An overview of SEACAS, including a short description of the codes in the system, will be presented. 
Abstracts and references for the codes are listed at the end of the report. 

Additional information about obtaining SEACAS can be obtained by contacting: 

Marilyn K. Smith 

Division 1425 

Sandia National Laboratories 

P.O. Box 5800 

Albuquerque, NM 87185-5800 

(505) 844-3082; Fax (505) 844-9297 


PR60IOING PAUF BLANK NOT FILMED 
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PHILOSOPHY 


SEACAS is a modular system based upon a common binary datafile format called EXODUS that 
includes the mesh description and the timeplanes of the computed results. A subset of this format, called 
GENESIS, is used to refer to the mesh description portion of the EXODUS format. 

All of the preprocessing, analysis, postprocessing, and translation codes can read and/or write 
EXODUS database files. A schematic of this is shown in Fig. 1. 


r 


Modular system - Each code tailored for a single function 

Two- and Three-Dimensional, Quasistatic and Dynamic Non- 
lear Analysis Codes 

Vectorized, Microtasked (Parallel execution 1-8 CPUS) 

/ 



Figure 1 - Modular structure of SEACAS 
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With this structure, codes can be tailored for a single function. For example, an analysis code can be 
added to the system without writing new mesh generation and postprocessing programs. Also, code-to- 
code data transfers and restarts use the EXODUS format. A more complete view of this is shown in Fig. 

2 . 


This modular concept is used extensively in the mesh generation process. Several special-purpose 
codes have been written that perform a certain mesh generation task. Examples of this are. 

• generate two-dimensional mesh (fastq) 

• transform two-dimensional mesh into a three-dimensional mesh (gen3d, genshell) 

• join two or more 2D or 3D meshes into a single mesh (gjoin) 

• reposition a 2D or 3D mesh (grepos). 

Finite element meshes of several complicated three-dimensional geometries have been successfully 
generated using this system, which is built up of codes that singly provide a limited capability, but when 
used as a system are very powerful. 

All of the codes are written in as portable a form as possible. Fortran codes are written in i ANSI 
Standard FORTRAN-77 and C language codes are written in ANSI Standard C where possible. Machine- 
specific routines are limited in number and are grouped together to minimize the time required to adapt 
them to a new system. 

A code management system is used for all of the SEA CAS codes to provide traceability and 
retrievability for quality assurance. The change logs include who changed the code, when it was chang , 
what was changed, and why the change was made. If required, a previous version can be retrieved. 



f BLOT 
FASTQ 
GEN3D 
GJOIN 
GREPOS 
Aprepro 
. PATRAN j 


GENESIS 


r PR0NT03D\ 
PRONTQ2D 
SANTOS2D 
SUBWAY 
SANTOS3D 
JAC2D 
JAC3D 
JACQ3D 
COYOTE 
COYOTE3D 
,ABAQUS j 
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Figure 2 - Ccxle-to-code data transfers using SEACAS 
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Figure 3 shows an example of a very large, complex problem solved using the codes in SEACAS. 

S EACA5 is divided into five broad categories of codes: prepost, analysis, translator, libraries, and 
scripts. The categories can be roughly defined as follows: 


prepost: Pre- and postprocessing codes including mesh generation, visualization, preprocessors 
and database manipulation codes. 


analysis: Finite element analysis codes including quasistatic, transient dynamics, thermal, and 
electromechanics. Two- and three-dimensional, nonlinear, large deformation. 


translators: Translation codes for editing output files (EXODUS), inter-machine translation, and 
exodus from/to commercial database translation. 


libraries: Support libraries including database routines, common machine-specific routines 
routines, graphics device drivers, and interactive help routines. 


plot 


scripts: Unix shell scripts for executing the prepost, analysis, and translation codes. Also includes 
support and installation routines. 


r 


425 fps Lateral Impact onto Rigid Surface 




240.000 eight-node hex elements 

290.000 Nodes (870,000 DOF) 

1 6 Materials 

19 Contact Surfaces 
PRONT03D Dynamic Analysis 
Microtasked on Cray Y-MP/864 
30 MWords memory 
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Figure 3 - Example of a complex problem solved using the codes in SEACAS 
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EXODUS DATA FILE 


Figure 4 summarizes the features of the EXODUS data file. An EXODUS data base contains data 
entities used for input and output of finite element analyses. The features of this data base are: 

• self-describing data file - the file includes information that describes the contained data 
(dimensionality, size, type, etc.). 

• random read (and write) access - data can be written or read in an arbitrary order. 

• machine-independent binary representation - data is stored in eXtemal Data Representation (XDR) 
format and thus can be transferred to, or accessed from, any machine without being concerned 
about what machine generated the data. 

• FORTRAN and C subroutine interface - application programs access the data via simple calls to a 
subroutine library. 

• extensible - because the file is self-describing and not just a rigid file format, new data entities and 
features can be added without modifying all of the application codes using the data base. 

• translators - capability to translate to/from many commercial application codes (i.e., ABAQUS, 
PATRAN, etc.). 




Common binary datafile format 

One common data file that each analysis code uses to communicate 
between the pre- and post-processors. 

• Random read (and write) access 

• Machine-independent binary representation 

• Fortran and C callable subroutine interface 

• Object-oriented data 

• Extensible 

•Translators available for PATRAN, ABAQUS, ASCII 

J 


V 


Figure 4 - Summary of Exodus data file 
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The EXODUS II data model is used for transferring finite element analysis data among application 
codes. It includes data to define the finite element mesh and label both boundary condition and load 
application points. EXODUS II accommodates multiple element types and is sufficiently general to service 
many different analysis codes, providing a very general format for analysis results. The data is stored in a 
machine-independent format and is randomly accessible. A FORTRAN and C application programming 
interface is also specified. 6 

Translation codes are used to convert EXODUS databases into different formats and to edit EXODUS 
databases. 

algebra. Manipulates EXODUS finite element output data by evaluating algebraic expressions. 
Equation variables are dependent on the input database variable names. 

seaexo: Converts SEACO files (the binary file format that preceded EXODUS) into EXODUS files. 

exotxt and txtexo: Converts EXODUS to and from a specially formatted ASCII text file. 

exoexo: A program which simply translates an EXODUS file to the same EXODUS file. It is used 
as a base program for writing new translators or database manipulation programs. 

abaexo: Converts ABAQUS (commercial finite element code) results files to EXODUS. 

con ex. Concatenates several EXODUS files into a single EXODUS file. Used to create a single 
EXODUS file from analyses that have been restarted. 8 

exoxdr and xdrexo: Converts EXODUS to and from an external data representation format (XDR) 
which can be transmitted between computers of different architectures, word lengths and byte 
orders. 

exogen: Creates a GENESIS mesh database from a specific time step of an EXODUS file. Used 
when an analysis of a deformed geometry is required. For example, an impact analysis followed by 
a thermal analysis. 

merlin2: Transfers (maps) nodal data between finite element meshes. For example, thermal output 
from a heat conduction code to a thermal input file for a quasistatic mechanics code. 

exopat and patexo: Converts EXODUS to and from PATRAN (commercial mesh generation 
program) neutral file format. 

exolexo2 and exo2exol: Converts EXODUS I format files to EXODUS II format. 

In addition to the translation codes listed above, the Cray Unicos system provides the capability to 
translate EXODUS files into a VMS and IEEE format using the exoexo code. 
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PREPROCESSING AND POSTPROCESSING CODES 


The pre- and postprocessing codes are comprised of the following mesh generation, visualization and 
database manipulation codes. 

gjoin: Join together two or more GENESIS databases into a single GENESIS database. 

gen3d: Transform a two-dimensional GENESIS database into a three-dimensional GENESIS 
database. Several transformations are supported and additional transformations can be easily 
added. 

grepos: Reposition or scale a GENESIS database. 

fastq: Interactive two-dimensional finite element mesh generation program. Includes several 
mesh generation options including paving. 

aprepro: An algebraic preprocessing program. 

genshell: Transform a two-dimensional GENESIS database into a three-dimensional shell 
GENESIS database. Several transformations are supported and additional transformations can be 
easily added. 

numbers: Calculates several properties of an EXODUS file, including mass properties, timesteps, 
condition numbers, cavity volumes, and others. 

grope: Interactively examine an EXODUS database. Grope is also contained within the blot 
program. Grope is primarily used to validate EXODUS files. 

blot: The primary graphical two-dimensional and three-dimensional postprocessing code. It includes 
deformed mesh plots, contour plots, shaded fringe plots, variable versus variable and time history 
plots, and distance versus variable plots. 


215 



Paving Technique 

Figure 7 shows examples of the paving technique. The paving algorithm fills an area by first placing 
elements along the interior and exterior boundaries. Row endpoints are determined so that elements can be 
fitted into and around comers. Once the row is generated, it is smoothed. The new row can be adjusted 
based on curvature. Elements called wedges are inserted into rows that are being generated concave 
outward to keep element sizes from becoming larger. Tucks or the removal of elements are performed for 
rows that are concave inward to keep elements from becoming smaller. When two rows of elements 
overlap, the overlap is sewn together by being pushed apart and attached together. When two rows are 
close but not touching, they are pulled together with a process called seaming cracks. The meshing 
process produces nearly square elements with comers of approximately 90 degrees. Smoothing and 
deletions of elements are performed until each element passes a "goodness" test. 


f The Paving Algorithm ' 

® 

Generate Boundaries 
Choose Row Endpoints 
Generate a Row 
Adjust the new Row 

• Wedges, Tucks and Smooth 

Seam Cracks 
Sew Overlaps 

\ J 



Figure 7 - Examples of the paving technique 
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PREPROCESSING AND POSTPROCESSING CODES 


The pre- and postprocessing codes are comprised of the following mesh generation, visualization and 
database manipulation codes. 
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grope: Interactively examine an EXODUS database. Grope is also contained within the blot 
program. Grope is primarily used to validate EXODUS files. 
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deformed mesh plots, contour plots, shaded fringe plots, variable versus variable and time history 
plots, and distance versus variable plots. 
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Figure 5 summarizes the capability the new pre-processing tool APREPRO gives a user. We have 
found that when APREPRO is combined with the other pre-processing tools, a great increase in 
productivity can be obtained. APREPRO can also be used to convert data in a variety of different units to 
one set of units. 


f " 'S 

APREPRO: Algebraic Preprocessor 

• Define problem in terms of variables 
•Global variables in one file that is included in all other files 
•Trigonometric and algebraic functions — C-like syntax 
•Portable: versions available on Unicos/VMS/Ultrix/Unix/Amiga 

Easily extendable (we wrote, have source) 


INPUT 


OUTPUT 


{ include ( one . inc ) } 

point 1 {xl = rad*cosd(angle) } {yl=rad} 
point 2 {x2 = xl + rad} {yl = yl + 1/2} 

{units (inch-lbf.-sec) } 
point 3 <5*m} {2.54*cm} 
point 4 {l*ft} {l*ft + l*in} 
velocity = {l*mile/hour} 


S Aprepro (version #) date 
§ include file ONE. INC 
$ rad = 1. 

$ angle = 30. 
point 1 0.866025 1.0 
point 2 1.866025 1.5 

$ units: inch- lbf -sec 
point 3 196.85039 1.0 
point 4 12.0 13.0 
velocity = 17.6 


• Any expression inside { } will be evaluated and printed 

• Can be used with any file that does not contain { } for other reasons 




Figure 5 - Aprepro: An algebraic preprocessor 
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MESH GENERATION 


Figure 6 lists the current areas of active research on mesh generation at Sandia. The r ^t goal of the 
CUBIT mesh generation project is to develop advanced meshing algorithms for meshing arbi ary 
JS£!i /voiumes. The pavfng algorism was develop a 2D 
algorithm It was recendy extended to meshing general 3D trimmed surfaces. The plastering algorithm is 
being developed as a 3D all-hexahedral volume meshing algorithm. 

The development of all of the 3D meshing algorithms cannot progress rapidly^ "ithouuh e inclusion of 
a general geometry generation mechanism. A solid modeling package, ACIS, is being used for t 
purpose. A nonmanifold topology and meshing database was designed and coded. This database serves 
all the mesh generation and geometry generation tools. The database allows direct access to 1 § 
functionality 8 and should prove an extremely useful prototype for commercial development of a true solid- 

model-based meshing program. 

The second goal is to develop for the first time an extremely robust adaptive finite element analysis 
capability for use in all fields of mechanics. The paving and plastering algorithms are very well suited for 

use in adaptive solution algorithms. 

The mesh generation algorithms for paving and plastering are superior to any mesh ^ ner ahon 
algorithms currently available to private industry. Consequently, seventi software houses have entered 
into Cooperative Research and Development Agreements (CRADA) with Sandia. 


CUBIT Mesh Generation Project 


• Algorithm Development 

• PAVING 

• PLASTERING 

• Solid Modeling 

• User Interface 

• Adaptive Analysis 

• Industry Consortium 



Figure 6 - Current areas of active research on mesh generation at Sandia 


217 




Paving Technique 


Figure 7 shows examples of the paving technique. The paving algorithm fills an area bv first niacin p 

filTnl a °H g 1,16 in , ten ° r and e o Xterio [ toundaries - Row end P°ints Je determine^soAatele^nte^^ 
fitted into and around comers. Once the row is generated, it is smoothed. The new row can be adjusted 
based on curvature. Elements called wedges are inserted into rows that are being generated concave 
outward to keep element sizes from becoming larger. Tucks or the removal of elements are performed for 
rows that are concave inward to keep elements from becoming smaller. When two rows of elements 
overlap, the overlap is sewn together by being pushed apart and attached together. When two rows are 
close but not touching, they are pulled together with a process called seaming cracks. The meshing 
process produces nearly square elements with comers of approximately 90 degrees. Smoothing and 
deletions of elements are performed until each element passes a "goodness" test. 




The Paving Algorithm 


—min. 1 1 


o 

— 

xmnxm 




v. 


Generate Boundaries 
Choose Row Endpoints 
Generate a Row 
Adjust the new Row 
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Figure 7 - Examples of the paving technique 
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As an example, one of the meshing tools is the sewing process. First intersections are found and then 
candidates for connection are chosen. Then the connections are formed between nodes that are close 
together. Seams are then formed by pushing rows of element edges together to remove the overlaps. 
Seams are then formed by pushing rows of element edges together to remove the overlaps. The mesh is 
then smoothed and elements not meeting a set of "goodness" tests are deleted. Figure 8 shows how 
sewing can be used to mesh a part with multiple holes. 
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Plastering Algorithm 


Figure 9 shows an example of a plastered volume. The plastering algorithm will be capable of filling 
an arbitrary three-dimensional solid with an all-hexahedral mesh. The algorithm will first mesh the 
exterior and interior surfaces with the surface paving algorithm. Then the volume will be filled by 
projecting the surface elements toward the interior to form elements. Once this is done to all the surfaces, 
this will define a new volume to be meshed. The process is continued until the volume is filled. Concepts 
analogous to sewing, seaming, wedges and tucks for surface paving will be used to place elements on the 

interior of surfaces. This is an example test case for a transition from a 2 x 2 mesh to a 4 x 4 mesh with 
nonregular meshing of the other sides. 


r 


Transition Model 


PAVED Surfaces PLASTERED Volume 





Figure 9 - Example of a plastered volume 
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Figure 10 shows 3D paving of a NURBS surface (nonuniform rational B-spline). The development of 
the surface paving algorithm has progressed to being capable of paving NURBS surfaces. This is an 
example of the surfaces being defined by a solid model and the paving algorithm meshing the surfaces. 

The example demonstrates the ability to pave arbitrary non-planar surfaces that have complex intersections 
with other surfaces. 



Figure 10 - Three-dimensional paving of a NURBS surface 
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Visualization Software 

Figure 15 shows a typical screen of what an analyst might expect to see while working at his desk. 
Application Visualization System (AVS), a commercial scientific visualization environment, was selected 
as the foundation for our software development effort after an extensive survey of the market two years 
ago. It was chosen because of the following characteristics: 

• functional -- contains a full-featured suite of modules which operate on numerous data types, 
including two-dimensional images and scalar and vector fields associated with structured or 
unstructured grids. 

• extensible - due to its modularity, functions can be readily added/customized. 

• ubiquitous — available on wide variety of hardware platforms. 

• distributable - AVS modules can execute on different machines within an application. 



Figure 15 - Example of the AVS visualization tool 
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Figure 10 shows 3D paving of a NURBS surface (nonuniform rational B-spline). The development of 
the surface paving algorithm has progressed to being capable of paving NURBS surfaces. This is an 
example of the surfaces being defined by a solid model and the paving algorithm meshing the surfaces. 

The example demonstrates the ability to pave arbitrary non-planar surfaces that have complex intersections 
with other surfaces. 



Figure 10 - Three-dimensional paving of a NURBS surface 
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Adaptive Analysis 


The paving algorithm has been used in a proof of concept example to perform an adaptive analysis. 
The geometry and loads for the problem are first defined. Then an initial mesh of approximately equal size 
elements is generated with the paving algorithm, and a static analysis is performed. The resulting stresses 
are used to determine an error function at the nodes of the initial mesh. The surface is then remeshed with 
the paving algorithm by sizing and placing elements based on an error function. This results in smaller 
elements where the error is large. This process was repeated four times and the final mesh produces an 
accurate answer to the problem. An example of adaptive meshing using paving is shown in Fig. 11. 



Figure 1 1 - Adaptive proof of concept example 
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Mesh Generation Consortium 

Sandia has created a Mesh Generation Consortium of computer-aided software companies and Sandia 
National Laboratories to develop and commercialize the meshing algorithms. To date, ten companies have 
shown active interest and Cooperative Research and Development Agreements (CRADAs) have been 
completed. The work is supported with industry funds and matching Department of Energy research 
funds. The paving algorithm will appear in commercial software within three months. Figure 12 lists the 
important points of the mesh generation consortium. 



Technology Transfer 
Industry Consortium 

Mesh Generation Consortium 

•Ten Companies with active interest 
•Four CRADA agreements are complete 
•Ford Motor Company 

• Fluid Dynamics International 
•MacNeal Schwendler Corporation 

• PDA Engineering 

• Matching DOE funds 

• PAVING to be in commercial software within three months 

• FORD has implemented PAVING in PDGS 



v. 



Figure 1 2 - Mesh generation consortium 
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PRODUCTION VISUALIZATION ENVIRONMENT 


Applied Visualization Group 

Supercomputers, both traditional and massively parallel machines, are capable of generating enormous 
amounts of data. Scientific visualization techniques are necessary to obtain useful information from these 
massive amounts of data. These techniques include responsive manipulation and animation of three- 
dimensional objects. 

Resulting databases can be on the order of several gigabytes, so high speed access to these large 
remote databases is required. Another challenge is that analysts want to perform these visualization 
activities at their desks. 

An Applied Visualization Group (AVG) has been established at Sandia National Laboratories to 
design and implement a production scientific visualization environment for use by the staff in the 
engineering sciences discipline. These analysts perform finite element and finite difference calculations in 
the areas of high velocity impact physics, shock wave physics, structural dynamics, structural mechanics, 
thermodynamics, and fluid mechanics. Visual results are critical in the data analysis functions performed. 


r 






Visualization System for Engineering Mechanics 


•A VS Software 
•Video Animation 
•Visualization server 
•Large Databases 


Super- 
computer Visualization 



•High Speed Communications 






Figure 13 - Visualization system for engineering mechanics 
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Visualization Server 

A visualization server concept is being implemented. One or several very high powered graphics 
machines will be used to perform the graphics manipulations, with the resulting images transmitted to an 
X-window display on the analyst's desk. We also expect to use PEX to distribute 3D graphics between 
the visualization server and the desktop display. We anticipate having a 100MByte channel connection 
from the visualization server to a supercomputer and network storage system, with FDDI connections 
from the visualization server to the offices. A block diagram of the proposed production environment 
consisting of hardware, software, and output capabilities integrated in an easy to use fashion is shown in 
Fig. 13. (See previous page.) 

A video animation system, using a component video laser disk recorder, high resolution converter, 
video transcoder, and a workstation, allows workers to rapidly create videos of analyses. The use of 
component and composite video technology allows for high quality when results are played directly from 
the laser disc, while still allowing for VHS video tapes to be created for showing at remote locations. A 
block diagram of the system currently in use at Sandia is shown in Fig. 14. 
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Figure 14 - Video animation system at Sandia 
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Visualization Software 

Figure 15 shows a typical screen of what an analyst might expect to see while working at his desk. 
Application Visualization System (AVS), a commercial scientific visualization environment, was selected 
as the foundation for our software development effort after an extensive survey of the market two years 
ago. It was chosen because of the following characteristics: 

• functional — contains a full-featured suite of modules which operate on numerous data types, 
including two-dimensional images and scalar and vector fields associated with structured or 
unstructured grids. 

• extensible — due to its modularity, functions can be readily added/customized. 

• ubiquitous — available on wide variety of hardware platforms. 

• distributable — AVS modules can execute on different machines within an application. 



Figure 15 - Example of the AVS visualization tool 
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ANALYSIS PROGRAMS 


The analysis codes in SEACAS include quasistatic, transient dynamics, thermal, and electromechamcs 
codes. 

iac2d: Jac2d is a finite element program which uses a nonlinear conjugate gradient technique to solve 
the large displacement, large strain, temperature-dependent and material nonlinear problem for two- 
dimensional plane or axisymmetric solids. 

iac3d: Jac3d is a finite element program which uses a nonlinear conjugate gradient technique to solve 
the large displacement, large strain, temperature-dependent and material nonlinear problem for three- 
dimensional solids. 

jacq3d: Jacq3d is a finite element program which uses a nonlinear conjugate gradient technique to 
solve the thermal conduction problem for three-dimensional solids. 

pronto2d: Pronto2d is a two-dimensional (planar and axisymmetric) transient solid dynamics code. 
Lagrangian formulation with explicit time integration is used for analyzing large deformations of 
highly nonlinear materials subjected to extremely high strain rates. 

pronto3d: Pronto3d is a three-dimensional transient solid dynamics code. Lagrangian formulation 
with explicit time integration is used for analyzing large deformations of highly nonlinear materials 
subjected to extremely high strain rates. 

sancho: Sancho is a finite element program which uses a dynamic relaxation technique to compute 
the quasistatic, large deformation, temperature-dependent, inelastic response of planar or axisymmetric 

solids. 

santos: Santos is a finite element program which uses a dynamic relaxation technique to compute the 
quasistatic, large deformation, temperature-dependent, inelastic response of planar or axisymmetric 
solids. Newer than sancho. 

santos3d: Santos3d is a three-dimensional version of santos. It is currently in development. 

conchas: Conchas is a linear structural analysis code for axisymmetric structures with loads that are 
symmetric about a plane. 

subway: Subway is a three-dimensional finite element program for calculating the transient electro- 
mechanical response of dielectric materials. 
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Constitutive Models 

Many constitutive models are implemented in the analysis codes. 

Elastic: Typical linear elastic material using Hooke's Law. 

Elastic-Plastic with Combined Hardening: Standard von Mises type yield condition with 
combined kinematic and isotropic linear hardening. 3V 3 

Viscoplastic: Simple rate-dependent plasticity. 

Viscoelastic: Thermoviscoelastic model for glass solidification. 

Damage: Dynamic fracture of brittle rock. 

Soils and Crushable Foam: Volumetric plasticity model. 

Low Density Foam: Low density polyurethane foam behavior. 

Hydrodynamic: Used with Equations of State in PRONTO 

Elastic-Plastic Hydrodynamic: Combination of elastic-plastic combined hardening with 
^uat^ofStote 6551116 response: Mie - Gruneisen type equation of state, JWL High Explosive 

Rate and Temperature Dependent Plasticity: Unified creep plasticity. 

puchWPi ry i :r D ep: Power harde ™ n g steady-state creep with elastic bulk response. Isotropic 
Elastic/Plastic Power- Law Hardening with Luders Strain. 

Johnson-Cook Strength: Rate and temperature-dependent plasticity. 

Elastic Plastic Power Law Hardening with Luder’s Strain: Describes post-yield strain 
hardening by a power law equation and includes a yield plateau. 

Hyperelastic: Stress based on the principal strains. 

bulk response ati ° n: P ° WCr hardening stead y- sta te creep deviatoric response with consolidation 

a l 1 °[. A . e constitutive models are implemented in all of the analysis codes. The quasistatic analysis 
consmutive models using temperature-dependent properties to facilitate 
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TRANSIENT DYNAMICS DEVELOPMENT 


A summary of the transient dynamics code PRONTO and the development objectives for PRONTO are 
shown in Fig. 16. 





PRONTO Description: 

A two- and three-dimensional transient solid dynamics code for ana- 
lyzing large deformations of highly non-linear materials subjected to 
high strain rates 

•Explicit mid-point time integration 
•One point element integration 
•2D - four node quad 

•3D - eight node hexagon, 3D - four node shell 
•Hourglass control 

• Adaptive time step control 
•Symmetric contact algorithm 

• Objective material coordinate system 


Code Development Objectives: 

•Well documented and consistent 

• Accurate numerical algorithms which minimize error 

• Dependable and aborts executions by providing diagnostics 
•Executes rapidly 

• Uses existing pre- and post- processing software 

• Easy to add new constitutive models 



Figure 1 6 - PRONTO description and development objectives 
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New Contact Detection Algorithm 

An increasingly important aspect of large-scale finite element simulations is the efficient and accurate 
determination of contact between deformable bodies. Generally there are several aspects of the numerical 
model that make contact detection difficult. This section reviews some currently used contact detection 
techniques and outlines some difficulties associated with these algorithms. Then, a new algorithm is 
proposed which circumvents these difficulties. 

The key points of the new contact detection algorithm are: i) uses the existing contact penalty/kinematic 
constraint method (including partitioned contact); ii) easily models self-contacting surfaces; iii) 
automatically defines all surfaces given the mesh connectivity; iv) models tearing and eroding surfaces; v) 
reduces user input; vi) uses a fast, memory-efficient global search to decide what slave nodes are in 
proximity to a master surface; and vii) does a detailed contact check using projected movements of both the 
master and slave to determine the magnitude and direction of slave node penetration of the master surface. 
See Fig. 17. 



Contact surfaces 


New contact detection algorithm: 

•Uses the existing contact penalty/kinematic constraint method 
•Partitioned contact (both surfaces can be master and slave) 

• Static and dynamic friction 

• Self contacting surfaces 
•Automatically defines contacting surfaces 

• Allows for tearing and eroding surfaces 

• Reduces user input 

•Does not require slideline pairing 

• Performs a global search for contact 

• Efficient contact search (knlogw) 

• Resolves contact corner ambiguities 


■\ 




Figure 17 - Contact surfaces 
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Multibody Impact: Elastic-Plastic Bar Impacting Bricks 

The following example demonstrates two features of the contact detection algorithm that are important 
to this particular class of problems. The example shown in Fig. 18 has an elastic-plastic bar impacting a 
stack of 17 elastic bricks. A stationary elastic-plastic wall is also resting against the stack of bricks. 

An important feature of the contact detection algorithm is the spatial independence of the kn Log n 
vectorized sorting algorithm. The speed of the algorithm is independent on the location of the bricks.. This 
becomes particularly important late in this example since the volume that the bodies occupy is increasing. 

Another feature of the contact detection algorithm that is both more efficient and useful is that the 
surfaces are not required to be paired. In fact, the surfaces were automatically defined using a new surface 
definition algorithm. For problems where random contact is anticipated, as in this example, each body 

could potentially impact any other body. For a contact-pairing algorithm, this implies that n 2 b contact pairs 
would be necessary, with each pair having 2m slave nodes. For the current global contact searching 
algorithm, it would imply one search with slave nodes. Consider in this example 19 2 pairs would 

have to be defined since random contact is expected. Assuming that each block has m ~ 50 slave nodes, 
19 2 pairs would require 19 2 (2m log (2m)) = 239,843 comparisons, whereas the current global contact 
searching algorithm would require only 19m/o£ 2 (19m) = 9397 comparisons. 


Contact Detection Algorithm: Multi-Body Impact 

«ik, 





1^' 

V 


J 


Figure 18 - Elastic-plastic bar impacting bricks 
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Self-Contacting Impact: Buckling of Shell-Like Structures 

The following example demonstrates the self-contacting capability of the contact detection algorithm. 
This feature is important to modeling crash dynamics where buckling, tearing and self-contact is common. 
The example shown in Fig. 19 has an elastic plate impacting an elastic-plastic can (shell-like structure). 
The can is 0.7 inches thick and has an inside radius of 12.5 inches. After 10 milliseconds the can is nearly 
crushed flat with numerous folds and buckles that self contact. 


PR0NT03D Contact Detection 
Can Crush 

® 


Oblique impact of steel can at 285 mph. 



Figure 19 - Self-contacting impact of a buckling shell- like structure 
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Automatic Contact Surface Redefinition: Tearing/Petalling of a Plate 

In the following example, the capability of automatically redefining the contact surface is essential. 
This feature is also important for modelling crash dynamics where buckling, tearing and self-contact is 
common. The example shown in Fig. 20 has a solid elastic sphere impacting a thin elastic-plastic plate. 
The plate is 0.25 inches thick and has an overall dimension of 12.0 inches by 14.0 inches. After a few 
milliseconds the plate is damaged such that two tears develop orthogonal to one another. The tears are 
actually a series of deleted elements where the damage has accumulated to 1.0. The newly created surface 
is automatically included in the contact algorithm by redefining the surface after any elements are deleted. 
At late time, the edges of the tears are in contact with the elastic sphere. 


PRONT03D 

Tearing - Sandia Damage 










Figure 20 - Automatic contact surface redefinition: tearing/petalling of a plate 
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Contact Algorithm Implementation: Reduced User Input 

Reduced user input is a convenience that is an outcome of the global searching and automatic surface 
definition capability. One can, for example, selectively include all surfaces of a body composed of a 
material and/or include only a subset of the surface. Figure 21 shows the typical input required for a 
problem. Reducing the input required from the user minimizes the number of mistakes, which minimizes 
cost. 



Contact Surface 


Reduced user input 


contact material 1 
contact material 2 
contact material 5 
contact surface 100 
contact data material 2, 


material 5, 


friction 


.1 



Surface 100, materials 1, 2 and 5 will be globally searched for con- 
tacts. Contacts between material 2 and 5 will have a coefficient of 
friction equal to 0. 1 



Figure 21 - Reduced user input 
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Contact Detection Algorithm: Capturing Slave Nodes for a Master Surface 

It is widely held that the location phase is the most time-consuming part of the contact detection 
algorithm. Obviously, the most robust approach would be to check every slave node against every master 
surface every time step. This exhaustive global searching approach requires nodal distance calculations on 
the order n z , where n is the number of nodes. Several algorithms have been proposed to speed up the 
location phase. 

The proposed contact detection algorithm uses a new global search algorithm. It uses 7n memory 
locations and requires ( k nLogn) comparisons and depends on only the number of entities (n) to be 
compared. It takes advantage of the known positions of the slave nodes and master surfaces as well as the 
predicted positions in the next time step. After the slave nodes are sorted by x, y and z coordinate, the 
master surfaces are processed sequentially. This processing involves bounding the space occupied by a 
master surface at its known location and its predicted location. Figure 22 shows a bounding box for a 
master surface over one time step. Clearly, a slave node could only make contact with this master surface 
if it is no further from the bounding box than the distance it could move in the time step. Therefore, any 
slave node inside the capture box should be considered for potential contact with the master surface. 



Figure 22 - Contact search bounding box 
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Contact Detection Algorithm: Slave Node/Master Surface Contact 

Unfortunately there are many ambiguous cases when determining contact between a slave node and a 
master surface. One such ambiguity arises because the surface normal is not continuous. This can result 
in either missed contacts or multiple solutions. In the proposed algorithm, the motion of the slave node 
and the master surface are considered and these ambiguities are avoided. As shown in Fig. 23, the time of 
contact and the point of contact is determined for each possible master surface that a slave node can 
penetrate. The algorithm solves for the time when a slave node and master surface all lie in the same 
plane. When the slave node can contact more than one master surface, the minimum time to contact 
determines the correct master surface contact. 

Further ambiguity is introduced when a contact algorithm performs detailed contact checks based only 
on the estimated (deformed) configuration. If an algorithm does not make use of information about how 
the slave node deforms to the estimated position, the algorithm must decide which contact violation to 
enforce. This is commonly done by determining which surface is most opposed to the slave node surface 
normal. When two surfaces are already in contact this so-called strength of contact check can be effective. 
However, when the surfaces are just coming into contact this type of static contact check cannot 
consistently predict the correct contact. 

In the proposed algorithm, the motion of the slave node and the master surface are considered and 
these ambiguities are avoided. In determining the push back direction, a distinction is made between 
convex and concave surfaces. The push back direction for a convex surface is determined simply by the 
minimum distance to the master surface. The push back direction can be along the master surface normal 
or it may be defined by the minimum distance to a vertex. For a concave surface the push back direction is 
always along the normal of the master surface that the slave node was previously in contact with. 





m 


m 
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FUTURE CODE DEVELOPMENT 


Figure 24 lists some of the active areas of code development and research at Sandia. 


r 



• Incorporate rezoning/remeshing algorithms 

• Constitutive model development 

• Orthotropic materials 

• Damage and tearing constitutive relations 

• Add Belytschko’s physical hourglass control 

•Add Rashid’s fully objective incremental formulation 

• Add rigid bodies, beams, springs and dash pots 

• Couple particle hydrodynamics with structural dynamics 

• C++ and Massively Parallel Computers 


v. 



Figure 24 - Future code development 
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uniform strain hexahedral elements are used in the finite element formulation. A number of new 
numerical algorithms which have been developed for the code are described in this report. An adaptive 
time step control algorithm is described which greatly improves stability as well as performance in 
plasticity problems. A robust hourglass control scheme which eliminates hourglass distortions without 
disturbing the finite element solution is included. All constitutive models in PRONTO are cast in an 
unrotated configuration defined using the rotation determined from the polar decomposition of the 
deformation gradient. An accurate incremental algorithm was developed to determine this rotation and 
is described in detail. A robust contact algorithm was developed which allows for the impact and 
interaction of deforming contact surfaces of quite general geometry. Numerical examples are presented 
to demonstrate the utility of these algorithms. 

L. M. Taylor and D. P. Flanagan, "PRONTO 3D, A Three-Dimensional Transient Solid Dynamics 
Program," SAND87-1912, Sandia National Laboratories, Albuquerque, New Mexico, March 1989. 

An external code interface is defined which allows other transient applications to communicate with the 
PRONTO family of finite element programs. This interface is written in ANSI FORTRAN and allows 
an independent author to specify requirements for an external code to PRONTO. The interface is 
written such that updates to PRONTO will not require modifications to the external code. 

L. M. Taylor and D. P. Flanagan, "An External Code Interface for the PRONTO Family of Transient 
Solid Dynamics Programs," SAND87-3003, Sandia National Laboratories, Albuquerque, New Mexico, 
September 1988 

PRONTO 2D and PRONTO 3D are two- and three-dimensional solid dynamics codes for analyzing 
large deformations of highly nonlinear materials subjected to high strain rates. This newsletter is 
issued to document changes to these codes. As of this writing, the latest version of PRONTO 2D is 
Version 4.5.6, and the latest version of PRONTO 3D is Version 4.5.6. 

This update of the two codes discusses two major modifications to the numerical formulations, three 
new constitutive models, and the additions and improvements of contact surfaces. Changes in file 
formats, other miscellaneous revisions, and the availability of PRONTO 2D and PRONTO 3D are also 
discussed. In addition, updated commands for PRONTO 2D are provided in Appendix A of this 
newsletter. 

S. W. Attaway, "Update of PRONTO 2D and PRONTO 3D Transient Solid Dynamics Program," 
SAND90-0102, Sandia National Laboratories, Albuquerque, New Mexico, November 1990. 

PRONTO 3D is a three-dimensional transient solid dynamics code for analyzing large deformations of 
highly nonlinear materials subjected to high strain rates. It is a Lagrangian finite element program with 
explicit integration of the equations of motion through time. This report documents the implementation 
of a four-node quadrilateral shell element into Version 6.0 of PRONTO 3D. 

This report describes the theory, implementation and use of a four-node shell element. Also described 
are the required architectural changes made to PRONTO 3D to allow multiple element types. Several 
test problems are documented for verification of the PRONTO 3D implementation and for 
demonstration of computational savings using shell elements for thin structures. These problems also 
serve as examples for the user. A complete, updated list of the PRONTO 3D input commands is also 
included. 

V. L. Bergmann, "Transient Dynamics Analysis of Plates and Shells with PRONTO 3D," SAND91- 
1182, Sandia National Laboratories, Albuquerque, New Mexico, September 1991. 

This update discusses modifications of PRONTO 3D tailored to the design of fast burst nuclear 
reactors. A thermoelastic constitutive model and spatially variant thermal history load were added for 
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FUTURE CODE DEVELOPMENT 


Figure 24 lists some of the active areas of code development and research at Sandia. 




• Incorporate rezoning/remeshing algorithms 

• Constitutive model development 

• Orthotropic materials 
•Damage and tearing constitutive relations 
•Add Belytschko’s physical hourglass control 

• Add Rashid’s fully objective incremental formulation 

• Add rigid bodies, beams, springs and dash pots 

• Couple particle hydrodynamics with structural dynamics 
•C++ and Massively Parallel Computers 




Figure 24 - Future code development 
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Couple Particle Hydrodynamics with Structural Dynamics 

Figure 25 shows an example problem using a combination of finite elements and particle elements. 

One major difficulty associated with the Lagrangian finite element method is modeling materials with no 
shear strength; for example, gases, fluids and explosives biproducts. Typically these materials can be 
modelled for only a short time with a Lagrangian finite element code. Tangling of the mesh will eventually 
lead to numerical difficulties such as negative element area or "bow tie" elements. Remeshing will allow 
the problem to continue for a short while, but the lack of shear strength causes instabilities that prevent a 
complete analysis. 


Smooth particle hydrodynamics is a gridless Lagrangian hydrodynamic technique. Requiring no 
mesh SPH can model material fracture, large shear flows, and penetration. SPH treats material as 
particles that have their masses smoothed in space. The density is computed at a point by summing the 
contributions of the smoothed particle masses in the vicinity of the point. SPH computes the strain rate 
and the stress divergence based on the nearest neighbors of a particle. The nearest neighbors are 
determined using an efficient particle sorting technique. 


Unique Capability to Couple Hydrodynamics and Structural Dynamics 

SPH - a gridless Lagrangian technique tor hydrodynamics 
PRONTO - a finite element transient dynamics code 
Example problem; Copper Block Impacting Water at 500 m/sec 



Figure 25 - Coupling particle hydrodynamics and structural dynamics 
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The SPH computational technique was embedded within the PRONTO computer code. SPH elements 
are modeled within PRONTO as elements that have only one node. Each element has the typical element 
variables associated with it (stress, strain, rotation, stretch, density, energy and other state variables). By 
using the existing PRONTO architecture, the SPH elements used the same constitutive equations as used 
by the quadrilateral elements in PRONTO. 

The embedding of the SPH method within PRONTO allows part of the problem to be modeled with 
quadrilateral elements while other parts of the problem are modeled with the gridless SPH method. The 
quadrilateral elements are coupled to the SPH elements through a contact-like algorithm. 
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Object-Oriented Code 

Figure 26 outlines the objectives for a new code that would combine solid, fluid and heat transfer 
problems into one code architecture. The goal of this effort is to apply object-oriented concepts to the 
design of a code architecture for solving problems in solid mechanics, fluid mechanics, and heat trans er. 
Consolidating these areas into a single package has the potential of reducing the time required for code 
maintenance and adding new features. Such an architecture should facilitate development of the capability 
to solve strongly coupled problems among these three areas of engineering science. The code will be 
implemented in the object-oriented language C++ to simplify code reuse and provide for extensibility. The 
code architecture will be designed to run efficiently on message passing MIMD computers and to take 
advantage of vector capabilities that exist on some machines in that class. 


^ Object-Oriented Finite Element Code 
Architecture for Massively Parallel 
Computers 

One code architecture for: 

• Solid mechanics problems 
•Fluid mechanics problems 
•Heat transfer problems 



C++ object-oriented architecture will: 

• simplify porting to new generations of supercomputers 

•ease implementation of new element types, solution algorithms 

•provide environment that encourages teamwork in software development 

•facilitate reuse of software 

•reduce effort for software maintenance 

•provide environment for solving large, strongly coupled problems 

— ' 


Figure 26 - Object-oriented finite element code architecture for massively parallel computers 
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ABSTRACTS AND REFERENCES 


References and complete abstracts for selected analysis, preprocessing, postprocessing, translation 
and library codes are listed in this section. The abstracts and references for translators are not listed. 


Algebra 

TTie Algebra program allows the user to manipulate data from a finite element analysis before it is 
plotted. The finite element output data is in the form of variable values (e.g., stress, strain, and 
velocity components) in an EXODUS database. The Algebra program evaluates user-supplied 
functions of the data and writes the results to an output EXODUS database which can be read by 
plot programs. J 


Amy P. Gilkey, "ALGEBRA - A Program that Algebraically Manipulates the Output of a Finite 
Element Analysis (EXODUS Version)," SAND88-1431, Sandia National Laboratories, Albuquerque 
New Mexico, August 1988. mm. 


Blot 

Blot is a graphics program for postprocessing of finite element analyses output in the EXODUS 
database format. It is command driven with free-format input and can drive any graphics device 
supported by the Sandia Virtual Device Interface. 

Blot produces mesh plots with various representations of the analysis output variables. The major 
mesh plot capabilities are deformed mesh plots, line contours, filled (painted) contours, vector plots 
of two/three variables (e.g., velocity vectors), and symbol plots of scalar variables (e.g., discrete 
cracks). Pathlines of analysis variables can also be drawn on the mesh. Blot's features include 
element selection by material, element birth and death, multiple views for combining several displays 
on each plot, symmetry mirroring, and node and element numbering. 

Blot can also produce X-Y curve plots of the analysis variables. Blot generates time-versus-variable 
plots or variable-versus-variable plots. It also generates distance-versus-variable plots at selected time 
steps where the distance is the accumulated distance between pairs of nodes or element centers, 

Amy P. Gilkey and John H. Glick, "BLOT-A Mesh and Curve Plot Program for the Output of a Finite 
1 Ana ysis ’ SAND88-1432, Sandia National Laboratories, Albuquerque, New Mexico, August 


i nnn c>‘ .UP^ ates to the Postprocessing Program BLOT," memo to distribution dated August 21 

1990, Sandia National Laboratories, Albuquerque, New Mexico. 


Conchas 

CONCHAS is a linear finite element structural analysis code which is specialized for axisymmetric 
structures. Loads and responses are limited to those symmetric about a plane which includes the 
symmetricaxis of the structure. CONCHAS will perform eigenanalysis, static analysis, and dynamic 
analysis. The element library includes consistent-mass shell, solid, and beam elements, nonlinear 

unh ties’ b Lu^^a^KA N FAS TO* ^ ost P rocess ' n ^ * s ava ^ a ble with the separately supported 
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William C. Mills-Curran and Dennis P. Flanagan, "CONCHAS Users Manual," SAND88-1006, 
Sandia National Laboratories, Albuquerque, New Mexico, June 1989. 


Exodus and Gene sis File Format 

William C Mills-Curran, Amy P. Gilkey and Dennis P. Flanagan, "EXODUS: A Finite Element File 
Format for Pre- and Postprocessing," SAND87-2977, Sandia National Laboratories, Albuquerque, New 
Mexico, September 1988. 

L. M. Taylor, D. P. Flanagan and W. C. Mills-Curran, "The GENESIS Finite Element Mesh File 
Format," SAND86-09 10, Sandia National Laboratories, Albuquerque, New Mexico, May 1986. 


Fasta 

The FASTQ code is an interactive two-dimensional finite element mesh generation program. It is 
designed to provide a powerful and efficient tool to both reduce the time required of an analyst to . 
generate a mesh, and to improve the capacity to generate good meshes in arbitrary geometries. It is 
based on a mapping technique and employs a set of higher-order primitives which have been 
developed for automatic meshing of commonly encountered shapes (i.e., the triangle, semi-circle, etc.) 
and conditions (i.e., mesh transitioning from coarse to fine mesh size). FASTQ has been designed to 
allow user flexibility and control. The user interface is built on a layered command level structure. 
Multiple utilities are provided for input, manipulation, and display of the geometric information, as 
well as for direct control, adjustment, and display of the generated mesh. Enhanced boundary flagging 
has been incorporated and multiple element types and output formats are supported. FASTQ includes 
the paving algorithm which meshes arbitrary two-dimensional geometries with an all quadrilateral 

mesh. 

T. D. Blacker, "FASTQ Users Manual, Version 2.1," SAND88-1326, Sandia National Laboratories, 
Albuquerque, New Mexico, July 1988. 

This paper presents a new mesh generation technique, paving, which meshes arbitrary two- 
dimensional geometries with an all-quadrilateral mesh. Paving allows varying element size 
distributions on the boundary as well as the interior of a region. The generated mesh is well formed 
(i.e., nearly square elements, elements perpendicular to boundaries, etc.) and geometrically pleasing 
(i.e., mesh contours tend to follow geometric contours of the boundary). In this paper we describe th< 
theory behind this algorithmic/heuristic technique, evaluate the performance of the approach and 
present examples of automatically generated meshes. 

T D. Blacker and M. B. Stephenson, "Paving: A New Approach to Automated Quadrilateral Mesh 
Generation," International Journal for Numerical Methods in Engineering, Vol. 32, 1991, pp. 811-847. 


Grope 

Grope is a program that examines the input to a finite element analysis (which is in the GENESIS 
database format) or the output from an analysis (in the EXODUS database format). Grope allows the - 
user to examine any value in the database. The display can be directed to the user s terminal or to a 
print file. 

Amy P. Gilkey, "GROPE - A GENESIS/EXODUS Database Examination Program," RS 1523/88/02, s 
Sandia National Laboratories, Albuquerque, New Mexico. * 


242 



Numbers 

Numbers is a program which reads and stores data from a finite element model described in the 
EXODUS database format. Within this program are several utility routines which generate information 
about the finite element model. The utilities currently implemented in Numbers allow the analyst to 
determine information such as (1) the volume and coordinate limits of each of the materials in t e 
model; (2) the mass properties of the model; (3) the minimum, maximum, and average element 
volumes for each material; (4) the volume and change in volume of a cavity; (5) the nodes or elements 
that arc within a specified distance from a user-defined point, line, or plane; (6) an estimate of the 
explicit central-difference timestep for each material; (7) the validity of contact surfaces or sidelines; 
that is, whether two surfaces overlap at any point; and (8) the distance between two surfaces. 

G. D. Siaardema, "NUMBERS: A Collection of Utilities for Pre- and Postprocessing Two- and 
Three-Dimensional EXODUS Finite Element Models," SAND88-0737, Sandia National Laboratories, 
Albuquerque, New Mexico, March 1989. 


Plasterine 

This report describes the progress of the three-dimensional mesh generation research, using plastering, 
during the 1990 fiscal year. Plastering is a three-dimensional extension of the two-dimensional paving 
technique. The objective is to fill an arbitrary volume with hexahedral elements. Tha plastering 
algorithm's approach to the problem is to remove rows of elements from the exterior o the vo ume. 
Elements are removed, one level at a time, until the volume vanishes. Special closure algorithms may 
be necessary at the center. The report also discusses the common development environment and 
software management issues. 

M. B. Stephenson, S. A. Canann and T. D. Blacker, "Plastering: A New Approach to Automated, 
Three-Dimensional Hexahedral Mesh Generation, Progress Report I," SAND89-2192, Sandia National 
Laboratories, Albuquerque, New Mexico, February 1992. 


Pronto2D and Pronto3D 

PR0NT02D is a two-dimensional transient solid dynamics code for analyzing large deformations of 
highly nonlinear materials subjected to extremely high strain rates. This Lagrangian finite element 
program uses an explicit time integration operator to integrate the equations of motion. Four node 
uniform strain quadrilateral elements are used in the finite element formulation. A number o new 
numerical algorithms which have been developed for the code are described m this report. An adaptive 
time step control algorithm is described which greatly improves stability as well as performance in 
plasticity problems. A robust hourglass control scheme which eliminates hourglass distortions without 
disturbing the finite element solution is included. All constitutive models in PRONTO are cast in an 
unrotated configuration defined using the rotation determined from the polar decomposition of the 
deformation gradient. An accurate incremental algorithm was developed to determine this rotation and 
is described in detail. A robust contact algorithm was developed which allows for the impact and 
interaction of deforming contact surfaces of quite general geometry. A number of numerical examples 
are presented to demonstrate the utility of these algorithms. 


Lee M. Taylor and Dennis P. Flanagan, "PR0NT02D, A Two-Dimensional Transient Solid Dynamics 
Program," SAND86-0594, Sandia National Laboratories, Albuquerque, New Mexico, March 198/. 

PR0NT03D is a three-dimensional transient solid dynamics code for analyzing large deformations of 
highly nonlinear materials subjected to extremely high strain rates. This Lagrangian finite element 
program uses an explicit time integration operator to integrate the equations of motion. Eight-node 
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uniform strain hexahedral elements are used in the finite element formulation. A number of new 
numerical algorithms which have been developed for the code are described in this report. An adaptive 
12 con jrol algonthm is described which greatly improves stability as well as perfoimance in ? 
p asncity problems. A robust hourglass control scheme which eliminates hourglass distortions without 
disturbing the finite element solution is included. All constitutive models in PRONTO are cast in an 
unrotated configuration defined using the rotation determined from the polar decomposition of the 

k^p^ShSi 1 a n u ccurate incre mental algorithm was developed to determine this rotation and 

is described in detail. A robust contact algonthm was developed which allows for the impact and 
interaction of deforming contact surfaces of quite general geometry. Numerical examples are presented 
to demonstrate the utility of these algorithms. p presented 

Pr„ c L n,m" T SAN r ™ n 7 d Pon & ',' PR P 1 ? T ? 3D ' A ^-Dimensional Transient Solid Dynamics 

gram, SAND87-1912, Sandia National Laboratones, Albuquerque, New Mexico, March 1989. 

PRONTTVfamllv *? defined which alI ™. s ? ther applications to communicate with the 

PRONTO family of finite element programs. This interface is written in ANSI FORTRAN and allows 

an independent author to specify requirements for an external code to PRONTO. The interface is 
written such that updates to PRONTO will not require modifications to the external code. 

Q Ta y lor and D. P. Flanagan, "An External Code Interface for the PRONTO Family of Transient 

SepfembeT^g'g ^' 115, SAND87 - 3003 ’ Sandia National Laboratories, Albuquerque, New Mexico 

h 1 ? 20 and PRONTO 3D are two- and three-dimensional solid dynamics codes for analyzing 
large deformations of highly nonlinear materials subjected to high strain rates This newsletter is g 
issued to document changes to these codes. As of this writing, lie latest version of PR0NT02D is 
Version 4.5.6, and the latest version of PRONTO 3D is Version 4.5.6. 

TTiis update of the two codes discusses two major modifications to the numerical formulations three 

d VC m °n e S ’ and the - addltlons and improvements of contact surfaces. Changes in’ file 
formats, other miscellaneous revisions, and the availability of PRONTO 2D and PRONTO 3D are also 

newsletter In add “ 10n ’ Updatcd commands for PRONTO 2D are provided in Appendix A of this 

lANTWrt nlm W Q y ’ K K - ° f P ?° k NT ° 2 P and PR0NT0 ® Transient Solid Dynamics Program," 
SAND90-0102, Sandia National Laboratones, Albuquerque, New Mexico, November 1990. 

PRONTO 3D is a three-dimensional transient solid dynamics code for analyzing large deformations of 

Sjf n °f e ^ mat f f l ah subjected t0 h *g h strain rates ; H is a Lagrangian finite element program with 
of a fr! 10 ^ ratlon , 0 f, the e 9 uatl ^ n ^ °I motion through time. This report documents the implementation 
of a four-node quadrilateral shell element into Version 6.0 of PRONTO 3D. V 

This report describes the theory, implementation and use of a four-node shell element. Also described 
are the required architectural changes made to PRONTO 3D to allow multiple element types Several 
test problems are documented for verification of the PRONTO 3D implementation and for 
emonstration of computational savings using shell elements for thin structures. These problems also 

Suded eXamP S 6 A C ° mp,ete> UpdatCd Hst ° f 016 PRONTO 3D input comSTs ako 


1 1 ®? r Smann, Transient Dynamics Analysis of Plates and Shells with PRONTO 3D " SAND91- 
1 182, Sandia National Laboratories, Albuquerque, New Mexico, September 1991. 

This update discusses modifications of PRONTO 3D tailored to the design of fast burst nuclear 
reactors. A thermoelastic constitutive model and spatially variant thermal history load were added for 
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this special application. Included are descriptions of the thermoelastic constitutive model and the 
thermal loading algorithm, two example problems used to benchmark the new capability, a user s 
guide and PRONTO 3D input files for the example programs. The results from PRONTO 3D 
thermoelastic finite element analysis are benchmarked against measured data and finite difference 
calculations. 

PRONTO 3D is a three-dimensional transient solid dynamics code for analyzing large deformations of 
highly nonlinear materials subjected to high strain rates. The code modifications are implemented in 
PRONTO 3D Version 5.3.3. 

D. S. Oscar, S. W. Attaway and J. D. Miller, "Modifications of the PRONTO 3D Finite Element 
Program Tailored to Fast Burst Nuclear Reactor Design," SAND9 1-0959, Sandia National Laboratories, 
Albuquerque, New Mexico, August 1991. 


SUBWAY 3D is a three-dimensional finite element code for numerical simulation of the transient 
electromechanical response of dielectric materials. The code uses a preconditioned conjugate gradient 
method to solve Poisson's equation governing the electric potential in the dielectrics. This field solver 
is embedded in a modified version of the finite element transient dynamics code PRONTO 3D and 
allows solution for the electric response at each time step used in the explicit integration of the 
equations of motion. The code is structured to allow flexibility in formulation of initial-boundary value 
problems by permitting specification of electrical conductors and allowing these conductors to be 
connected to electrical circuits isolated from the mechanical deformations. An algorithm for solution of 
these initial-boundary value problems is incorporated into the code with special material models to 
represent the response of ordinary dielectrics and electromechanically active dielectrics like 
piezoelectric and ferroelectric ceramics. The code has a wide range of applications and, in particular, 
can be used to calculate responses of shock activated power supplies, impact gages, and the change in 
electrical capacitance due to deformations, 

S. T. Montgomery - documentation not yet available. 


Sancho 

SANCHO is a finite element computer program designed to compute the quasistatic, large 
deformation, inelastic response of planar or axisymmetric solids. Finite strain constitutive 
theories for plasticity, volumetric plasticity, and metallic creep behavior are included. A constant 
bulk strain, bilinear displacement isoparametric finite element is employed for the spatial discretization. 
The solution strategy used to generate the sequence of equilibrium solutions is a self-adaptive dynamic 
relaxation scheme which is based on explicit central difference pseudo-time integration and artificial 
damping. A master-slave algorithm for sliding interfaces is also implemented. A theoretical 
development of the appropriate governing equations and a description of the numerical algorithms are 
presented along with a user's guide which includes several sample problems and their solution. 

Charles M. Stone, Raymond D. Krieg and Zelma E. Beisinger, "SANCHO, A Finite Element 
Computer Program for the Quasistatic, Large Deformation, Inelastic Response of Two-Dimensional 
Solids," SAND84-2618, Sandia National Laboratories, Albuquerque, New Mexico, April 1985. 
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Santos and Santos3D 


SANTOS is a finite element computer program designed to compute the quasistatic, large deformation 
inelastic response of planar or axisymmetric solids. SANTOS is based on the dynamics program 
PR0NT02D by L. M. Taylor and D. P. Flanagan. SANTOS utilizes a self-adaptive dynamic 
relaxation algorithm to achieve a quasistatic solution. The efficiency, speed, through vectorization, 
and state-of-the-art algorithms that Taylor and Flanagan built into PRONTO2D are maintained in 
SANTOS. The architecture of the code as well as the user interface is similar for both codes which 
improves code reliability and encourages the use of both codes by analysts. By utilizing the same 
material interface, the same constitutive models may be utilized by both codes which allows for 
coupling of the two codes. 

C. M. Stone - documentation not yet available. 


Jac2D. .Iac3D and .Tacq3D 

The nonlinear conjugate gradient procedure is employed in the computer program JAC2D to solve 
quasi-static nonlinear mechanics problems. A set of continuum equations which describe accurately 
nonlinear mechanics involving large rotation and strain are very conveniently used with the conjugate 
gradient method to solve the nonlinear problem. The method is exploited in a two-dimensional setting 
while using various methods for accelerating convergence. Sliding interface conditions are also 
implemented. A four-node Lagrangian uniform strain element is used with hourglass stiffness to 
control the zero energy nodes. Materials which can be modeled are elastic and isothermal elastic- 
plastic with combined kinematic and isotropic hardening. The program is vectorized to perform very 
efficiently on CRAY computers. Sample problems described are the bending of a thin beam, the 
rotation of a unit cube, and a pressurized and thermally-loaded sphere and cylinder. 

J. H. Biffle, "A Two-Dimensional Finite Element Computer Program for the Nonlinear Quasistatic 
Response of Solids with the Conjugate Gradient Method," SAND81-0998, Sandia National Laboratories 
Albuquerque, New Mexico, April 1984. 

JAC3D is a three-dimensional finite element program designed to solve quasi-static nonlinear 
mechanics problems. A set of continuum equations, which describe nonlinear mechanics involving 
large rotation and strain, is used. A nonlinear conjugate gradient method is used to solve the nonlinear 
equations. The method is implemented in a three-dimensional setting with various methods for 
accelerating convergence. Sliding interface conditions are also implemented. An eight-node 
Lagrangian uniform strain element is used with hourglass stiffness to control the zero energy modes. 
This release of the code contains elastic and isothermal elastic-plastic material models. Other material 
models can be added relatively easily. The program is vectorized to perform very efficiently on CRAY 
computers. Sample problems described are the bending of a thin beam, the rotation of a unit cube, and 
a pressurized and thermally loaded sphere. 


J. H. Biffle, "JAC3D - A Three-Dimensional Finite Element Computer Program for the Nonlinear 
Quasistatic Response of Solids with the Conjugate Gradient Method," in preparation. 

The nonlinear conjugate gradient procedure is employed in the computer program JAC3D to solve the 
steady-state and transient nonlinear heat conduction problem for solids in three dimensions. The 
program can also be used for other types of diffusion problems. The problem is formulated with the 
finite element method and employs an eight-node uniform gradient element with hourglass stiffness to 
control the zero energy modes. JACQ3D is highly vectorized to perform efficiently on the CRAY 
computer. Sample problems are included to verify the code and to provide examples of the use of the 
code. 
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J. H. Biffle, "JAC3D - A Three-Dimensional Finite Element Computer Program for Nonlinear Heat 
Conduction Problems with the Conjugate Gradient Method," in preparation. 


Merlin II 

The MERLIN II program is designed to transfer data between finite element meshes of arbitrary 
geometry. The program is structured to accurately interpolate previously computed solutions onto a 
given mesh and format the resulting data for immediate use in another analysis program. Data from 
either two-dimensional or three-dimensional meshes may be considered. The theoretical basis and 
computational algorithms used in the program are described and complete user instructions are 
presented. Several example problems are included to demonstrate program usage. 

D. K. Gartling, "MERLIN II - A Computer Program to Transfer Solution Data Between Finite 
Element Meshes," SAND89-2989, Sandia National Laboratories, Albuquerque, New Mexico, July 1991. 


Supes Library 

The Software Utilities Package for the Engineering Sciences (SUPES) is a collection of subprograms 
which perform frequently used non-numerical services for the engineering applications programmer. 
The three functional categories of SUPES are: (1) input command parsing, (2) dynamic memory 
management, and (3) system dependent utilities. The subprograms in categories one and two are 
written in standard FORTRAN-77, while the subprograms in category three are written to provide a 
standardized FORTRAN interface to several system dependent features. 

J. H. Red-Horse, W. C. Mills-Curran and D. P. Flanagan, "SUPES Version 2.1, A Software Utilities 
Package for the Engineering Sciences," SAND90-0247, Sandia National Laboratories, Albuquerque, New 
Mexico, May 1990. 


Gen3D 

GEN3D is a three-dimensional mesh generation program. The three-dimensional mesh is generated by 
mapping a two-dimensional mesh into three-dimensions according to one of four types of 
transformations: translating, rotating, mapping onto a spherical surface, and mapping onto a 
cylindrical surface. The generated three-dimensional mesh can then be reoriented by offsetting, 
reflecting about an axis, and revolving about an axis. GEN3D can be used to mesh geometries that are 
axisymmetric or planar, but, due to three-dimensional loading or boundary conditions, require a three- 
dimensional finite element mesh and analysis. More importandy, it can be used to mesh complex 
three-dimensional geometries composed of several sections when the sections can be defined in terms 
of transformations of two-dimensional geometries. The code GJOIN is then used to join the separate 
sections into a single body. GEN3D reads and writes two-dimensional and three-dimensional mesh 
databases in the GENESIS database format; therefore, it is compatible with the preprocessing, 
postprocessing, and analysis codes used by the Engineering Analysis Department at Sandia National 
Laboratories, Albuquerque, New Mexico. 

Amy P. Gilkey and Gregory D. Sjaardema, "GEN3D: A GENESIS Database 2D to 3D 
Transformation Program, " SAND89-0485, Sandia NationalXaboratories, Albuquerque, New Mexico, 
March 1989. 

This memo describes the changes that have been made to the GEN3D program since the manual 
(SAND89-0485) was published. The changes include: (1) new mesh generation options: spline, 
project, twist, interval, and rotcen transformations; (2) new mesh modification options: change 










J. H. Biffle, "JAC3D - A Three-Dimensional Finite Element Computer Program for Nonlinear Heat 
Conduction Problems with the Conjugate Gradient Method," in preparation. 


Merlin II 

The MERLIN II program is designed to transfer data between finite element meshes of arbitrary 
geometry. The program is structured to accurately interpolate previously computed solutions onto a 
given mesh and format the resulting data for immediate use in another analysis program. Data from 
either two-dimensional or three-dimensional meshes may be considered. TTie theoretical basis and 
computational algorithms used in the program are described and complete user instructions are 
presented. Several example problems are included to demonstrate program usage. 

D. K. Gartling, "MERLIN II - A Computer Program to Transfer Solution Data Between Finite 
Element Meshes," SAND89-2989, Sandia National Laboratories, Albuquerque, New Mexico, July 1991. 


Supes Library 

The Software Utilities Package for the Engineering Sciences (SUPES) is a collection of subprograms 
which perform frequently used non-numerical services for the engineering applications programmer. 
The three functional categories of SUPES are: (1) input command parsing, (2) dynamic memory 
management, and (3) system dependent utilities. The subprograms in categories one and two are 
written in standard FORTRAN-77, while the subprograms in category three are written to provide a 
standardized FORTRAN interface to several system dependent features. 

J. H. Red-Horse, W. C. Mills-Curran and D. P. Flanagan, "SUPES Version 2.1, A Software Utilities 
Package for the Engineering Sciences," SAND90-0247, Sandia National Laboratories, Albuquerque, New 
Mexico, May 1990. 


Geti3D 

GEN3D is a three-dimensional mesh generation program. The three-dimensional mesh is generated by 
mapping a two-dimensional mesh into three-dimensions according to one of four types of 
transformations: translating, rotating, mapping onto a spherical surface, and mapping onto a 
cylindrical surface. The generated three-dimensional mesh can then be reoriented by offsetting, 
reflecting about an axis, and revolving about an axis. GEN3D can be used to mesh geometries that are 
axisymmetric or planar, but, due to three-dimensional loading or boundary conditions, require a three- 
dimensional finite element mesh and analysis. More importantly, it can be used to mesh complex 
three-dimensional geometries composed of several sections when the sections can be defined in terms 
of transformations of two-dimensional geometries. The code GJOIN is then used to join the separate 
sections into a single body. GEN3D reads and writes two-dimensional and three-dimensional mesh 
databases in the GENESIS database format; therefore, it is compatible with the preprocessing, 
postprocessing, and analysis codes used by the Engineering Analysis Department at Sandia National 
Laboratories, Albuquerque, New Mexico. 

Amy P. Gilkey and Gregory D. Sjaardema, "GEN3D: A GENESIS Database 2D to 3D 
Transformation Program," SAND89-0485, Sandia National Laboratories, Albuquerque, New Mexico, 
March 1989. 


This memo describes the changes that have been made to the GEN3D program since the manual 
(SAND89-0485) was published. The changes include: (1) new mesh generation options: spline, 
project, twist, interval, and rotcen transformations; (2) new mesh modification options: change 
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material, change sideset, and change nodeset; (3) new mesh orientation option: scale; and (4) 
miscellaneous changes. 

G. D. Sjaardema, "Updates to the mesh generation program GEN3D," memo to distribution dated 
April 11, 1990, Sandia National Laboratories, Albuquerque, New Mexico. 


Giflin 

Gjoin is a computer program which is used to merge two or more GENESIS databases into a single 
GENESIS database. 

G. D. Sjaardema, "GJOIN: A Program for Merging Two or More GENESIS Databases," (memo to 
distribution dated June 19, 1991), Sandia National Laboratories, Albuquerque, New Mexico. 


Grepos 

GREPOS is a mesh utility program that repositions or modifies the configuration of a 2D or 3D mesh. 
Grepos can be used to change the orientation and size of a 2D or 3D mesh; change the material block, 
nodeset, and sideset IDs; or "explode" the mesh to facilitate viewing of the various parts of the model. 
Grepos also updates the EXODUS QA and information records to help track the codes and files used 
to generate the mesh. GREPOS reads and writes 2D and 3D mesh databases in the GENESIS 
database format; therefore, it is compatible with the preprocessing, postprocessing, and analysis codes 
in SEACAS. 

G. D. Sjaardema, "GREPOS: A GENESIS Database Repositioning Program," SAND90-0566, 
Sandia National Laboratories, Albuquerque, New Mexico, April 1990. 


Aprepro 

APREPRO is a translator that reads a text file containing both general text and algebraic expressions. 

It echoes the general text to the output file, along with the results of the algebraic expressions. The 
syntax used in APREPRO is such that all expressions between the delimiters '{' and are evaluated 
and all other text is simply echoed to the output file. 

G. D. Sjaardema, "Aprepro: An Algebraic Preprocessor for Input Files," memo to distributed dated 
May 18, 1990, Sandia National Laboratories, Albuquerque, New Mexico. 

The initial implementation of a units conversion capability in Aprepro has been completed. The units 
conversion is a very useful capability for analysts and it is required for the implementation of the 
MATS material database system. 

G. D. Sjaardema, "Implementation of Units Conversion in Aprepro," memo to distribution dated April 
24, 1992, Sandia National Laboratories, Albuquerque, New Mexico. 


Constituti ve Models 


B. J. Thome, "A Damage Model for Rock Fragmentation and Comparison of Calculations with 
Blasting Experiments in Granite," SAND90-1389, Sandia National Laboratories, Albuquerque, New 
Mexico, October 1990. 
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E. P. Chen, "A Computational Model for Jointed Media with Orthogonal Sets of Joints, SAND86- 
1122, Sandia National Laboratories, Albuquerque, New Mexico, March 1987. 

M. K. Neilsen, H. S. Morgan, R. D. Krieg, "A Phenomenological Constitutive Model for Low 
Density Polyurethane Foams," SAND86-2927, Sandia National Laboratories, Albuquerque, New Mexico, 
April 1987. 

R. S. Chambers, "A Viscoelastic Material Model for Computing Stresses in Glass," SAND90-0645, 
Sandia National Laboratories, Albuquerque, New Mexico, July 1990. 

G. D. Sjaardema and R. D. Krieg, "A Constitutive Model for the Consolidation of WIPP Crushed Salt 
and Its Use in Analyses of Backfilled Shaft and Drift Configurations," SAND87-1977, Sandia National 
Laboratories, Albuquerque, New Mexico, October 1987. 

C. M. Stone, G. W. Wellman and R. D. Krieg, "A Vectorized Elastic/Plastic Power Law Hardening 
Material Model Including Luders Strain," SAND90-0153, Sandia National Laboratories, Albuquerque, 
New Mexico, March 1990. 

J R. Weatherby, R. D. Krieg and C. M. Stone, "Incorporation of Surface Tension into the Structural 
Finite Element Code SANCHO," SAND89-0509, Sandia National Laboratories, Albuquerque, New 
Mexico, March 1989. 

S. W. Attaway, "A Local Isotropic/Global Orthotropic Finite Element Technique for Modeling the 
Crush of Wood," SAND88-1449, Sandia National Laboratories, Albuquerque, New Mexico, September 
1988. 

R. D. Krieg, "A Simple Constitutive Description for Soils and Crushable Foams," SC-DR-72-0883, 
Sandia National Laboratories, Albuquerque, New Mexico, 1987. 

D. J. Bammann, "An Internal Variable Model of Viscoplasticity," International Journal of Engineering 
Science, Vol. 22, 1984, pp. 1041-1053. 
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LIST OF BOOKS, MONOGRAPHS, SPECIAL ISSUES OF JOURNALS, AL& 
SURVEY PAPERS AND REPORTS ON 
COMPUTATIONAL TECHNIQUES FOR SIMULATING CRASH 

Howard S. Levine 

Weidlinger Associates, Los Altos, CA 
and 

Ahmed K. Noor 

University of Virginia, Hampton, VA 

Several books, monographs, conference proceedings and research reports have been pub- 
lished on the subject of crash dynamics. For the benefit of the readers of the proceedings, all the 
literature known to the compilers is listed subsequendy. The literature is divided into five cate- 
gories: books, monographs and special issues of journals; survey papers; reports; research papers; 
and commercial software systems. The references within each category are listed alphabetically. 

I. Books. Monographs and Special Issues of Journals 

1 . Davies, G. A. O., ed.: Proceedings of the International Conference on Structural Impact and 
Crashworthiness, vol. 1, Keynote Lectures, Elsevier, London, 1984. 

2. International Symposium on Structural Crashworthiness and Failure. Special issue of 
International Journal of Mechanical Sciences, Vol. 35, Nos. 3/4, 1993. 

3. Johnson, W.; and Mamalis, A. G., eds.: Crashworthiness of Vehicles. Mechanical 
Engineering Publications, Ltd., London, 1978. 

4. Jones, N.; and Wierzbicki, T., eds ..Structural Crashworthiness. Butterworths, London, 1983. 

5. Jones, N., ed.: Structural Failure Symposium, June 6-8, 1988. Int. J. Impact Eng., special 
issue, vol. 7, no. 2, 1988. 

6. Jones, N., ed.: Symposium on Structural Crashworthiness and Failure. Special issue of 
International Journal of Impact Engineering, Vol. 13, No. 2, May 1993. 

7 . Jones, N.; and Wierzbicki, T., eds.: Structural Crashworthiness and Failure. Elsevier Science 
Publishers, London, 1993. 


8. Khalil, T. B.; Ne, C. M.; Mahmoad, H. F.; and King, A. I., eds.: Crashworthiness and 
Occupant Protection in Transportation Systems - 1991. AMD Vol. 126/BED Vol. 19, ASME 
Winter Annual Meeting, Atlanta, GA, 1991. 

9. King, A.; and Khalil, T., eds.: Crashworthiness and Occupant Protection. AMD Vol. 106, 
BED Vol. 13, American Society of Mechanical Engineers, NY, 1989. 

10. Macauley, M.: Introduction to Impact Engineering. Chapman & Hall, London, 1987. 

11. Morton, J., ed.: Proceedings of the International Conference on Structural Impact and 
Crashworthiness, vol. 2, Conference Papers, Elsevier, London, 1984. 

12. Reid, S. R., ed.: Structural Failure Symposium, June 6-8, 1988, Int. J. Mech. Sci., special 
issue, vol. 30, nos. 3/4, 1988. 


13. Sabir, A. B.; and Niku-Lari, A., eds.: CRASH-93: Numerical Methods in Crash Simulation. 

Il lT International, Goumay-sur-Mame, France. 
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14. Saczalski, K.: Aircraft Crashworthiness. Univ. Press of Virginia, Charlottesville, VA, 1975. 

15. Schwer, L. E.; Salamon, L. E.; and Liu, W., eds.: Computational Techniques for Contact, 

Impact, Penetration and Perforation of Solids. AMD Vol. 103, ASME Winter Annual 
Meeting, San Francisco, CA, Dec. 10-15, 1989, American Society of Mechanical Engineers, 
NY, 1989. 6 

16. Tong, P., Ni, C. M., King, A. and Lantz, S., eds.: Symposium on Vehicle Crashworthiness 
Including Impact Biomechanics. AMD Vol. 79, BED Vol. 1, American Society of Mechanical 
Engineers, NY, 1986. 

17. Wierzbicki, T.; and Jones, N., eds.: Structural Failures. John Wiley & Sons, NY, 1989. 

II. Survey Papers 

1 . Anderson, C. E.; and Bodner, S. R.: The Status of Ballistic Impact Modeling. Int. J. Impact 
Eng., vol. 7, 1988, pp. 9-35. 

2. Belytschko, T.: On Computational Methods for Crashworthiness. Comput. & Struct , vol 
42, no. 2, 1992, pp. 271-279. 

3. Campbell, G. S.; and Lahey, R. T. C.: A Survey of Serious Aircraft Accidents Involving 
Fatigue Fracture, vol. 2. Rotary-Wing Aircraft. National Aeronautical Establishment, Ottawa, 
Ontario, NRC 21277, April 1983. (See also vol. 1, AD-A137-254.) 

4. Campbell, G. S. and Lahey, R. T. C.: A Survey of Serious Aircraft Accidents Involving 
Fatigue Fracture. Vol. I. Fixed-Wing Aircraft. National Aeronautical Establishment, Ottawa, 
Ontario, NRC 21276, April 1983. (See also vol. 2, AD-A137-255.) 

5. Carden, H. D.: Impact Dynamics Research on Composite Transport Structures. NASATM- 
83691, March 1985. 

6. Dallard, P. R. B.; and Miles, J. C.: Design Tools for Impact Engineers. Proceedings of the 
International Conference on Structural Impact and Crashworthiness. G. A. O. Davies, ed., 
vol. 1, 1984, pp. 369-382. 

7 . Farley, G. L.; Boitnott, R. L.; and Carden, H. D.: Helicopter Crashworthiness Research 
Program. NASA CP-2495, 1987. 

8. Fasanella, E. L.; Carden, H. D.; Boitnott, R. L.; and Hayduk, R. J.: A Review of the 
Analytical Simulation of Aircraft Crash Dynamics. NASA TM- 102595, Jan. 1990. 

9. Hayduk, R. J.; Thomson, R. G.; and Carden, H. D.: NASA/FAA General Aviation Crash 
Dynamics Program - An Update. Forum , vol. 12, no. 3, 1979, pp. 147-156. 

10. Hayduk, R. J.; Carden, H. D.; Fasanella, E. L.; and Boitnott, R. L.: Status of Analytical 
Simulation of Aircraft Crash Dynamics. 66th AGARD Structures and Materials Panel Meeting, 
Luxembourg, May 1988. 

1 1. Huculak, P.: A Review of Research and Development in Crashworthiness of General Aviation 
Aircraft: Seats, Restraints and Floor Structures. Aeronautical Note NAE-AN-64. NRC No. 
31334, Ottawa, National Research Council of Canada, Feb. 1990. 
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12. Hui, D.; and Jones, N.: Recent Advances in Impact Dynamics of Engineering Structures - 
1989. Winter Annual Meeting, American Society of Mechanical Engineers, Dec. 10-15, 1989, 
San Francisco, CA, AMD-vol. 105, AD-vol. 17, ASME, NY, 1989. 

13. Jarzab, W.; and Schwarz, R.: Crashworthiness of Aircraft Structures (U). Advisory Group 
for Aerospace Research and Development, Neuilly-Sur-Seine, France. Conference 
Proceedings of Energy Absorption of Aircraft Structures as an Aspect of Crashworthiness 
(66th), Luxembourg, May 1-6, 1988. Report AD- A2 12-606, Dec. 1988, pp. 17.1-17.13. 

14. McComb, Jr., H.; Thomson, R. G.; and Hayduk, R. J.: Structural Dynamics Research in a 
Full-Scale Transport Aircraft Crash Test. J. Aircraft, vol. 24, no. 7, July 1987. 

15. Poon, G: A Review of Crashworthiness of Composite Aircraft Structures. Aeronautical Note 
NAE-AN-63. NRC No. 31276, Ottawa, National Research Council of Canada, Feb. 1990. 

16. Saczalski, K. J.; Singley III, G. T.; Pilkey, W. D.; and Huston, R. L.: Aircraft 
Crashworthiness. University Press of Virginia, Charlottesville, VA, 1975. 

17. Thomson, R. G.; and Goetz, R. C.: NASA/FAA General Aviation Crash Dynamics Program 
- A Status Report. J. Aircraft, vol. 17, no. 8, Aug. 1980, p. 584. 

18. Thomson, R. G.; Carden, H. D.; and Hayduk, R. J.: Survey of NASA Research on Crash 
Dynamics. NASA TP-2298, April 1984. 

19. Thomson, R. G.; Carden, H. D.; and Hayduk, R. J.: Research at NASA on Crash Dynamics , 
Structural Impact and Crashworthiness, vol. 1, Keynote Lectures, G. A. O. Davies, ed., 
Elsevier Applied Science Publishers, London, July 1984, pp. 1-43. 

20. Thornton, P. M.; Mahmood, H. F.; and Magee, C. L.: Energy Absorption by Structural 
Collapse. Proc. of the International Conference on Structural Impact and Crashworthiness, G. 
A. O. Davies, ed., vol. 1, 1984, pp. 96-1 17. 

21. Widmayer, E.: A Structural Survey of Classes of Vehicles for Crashworthiness. Boeing 
Vertol Co. Report No. FRA/ORD-79- 1 3, U.S. Dept, of Transportation, 1979. 

22. Wittlin, G.; and Gamon, M. A.: A Literature Survey of Airborne Vehicles Impacting with 
Water and Soil; Head Injury Criteria and Severity Index Development of Computer Program 
KRASH. DOT/FAA/CT- 90/24, July 1992. 

III. Reports 

1 . A Study of U.S. Air Carrier Accidents 1964-1969. NTSB Report No. PB-21 1-054, 1972. 

2. Advisory Group for Aerospace Research and Development, Neuilly-Sur-Seine, France. 
Conference Proceedings of Energy Absorption of Aircraft Structures as an Aspect 
Crashworthiness (66th) held in Luxembourg, May 1-6, 1988 (U). Report AGARD-CP-443. 

3. Advisory Group for Aerospace Research and Development, Neuilly-Sur-Seine, France. 
Crashworthiness of Airframes (U). Presented at the 61st Meeting of the Structures and 
Materials Panel of AGARD, Oberammergau, Germany, Sept. 8-13, 1985. 

4. Ahlers, R. H.: Full-Scale Aircraft Crash Tests of Modified Jet Fuel. U.S. Dept, of 
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